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SUMMARY 

MODIFLY  is  a modular  trajectory  simulation  computer  program  v^ich  was  written 
so  as  to  be  effective,  efficient,  and  easily  modified.  The  program  was  designed 
primarily  for  the  simulation  of  typical  autonomous  guided  missiles  in  which 
roughly  equal  consideration  is  given  to  the  simulation  of  the  seeker,  guidance, 
autopilot,  controls,  and  aerodynaunics.  It  was  written  in  FORTRAN  IV  language 
for  use  on  a CDC  6500  computer  emd  uses  overlay  files  and  a library  edit  routine. 
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INTRODUCTION 

This  program  was  written  in  order  to  accomplish  two  specific  tasks.  The 
first  was  to  decrease  costs  and  the  second  was  to  increase  efficiencies  in  the 
simulation  of  the  flight  of  any  vehicle  that  moves  above  the  earth's  surface. 
These  goals  were  attained  by  preparing  a program  which  allowed  for  the  condensing 
of  several  older  NSWC  trajectory  programs  into  as  small  a package  as  possible 
in  order  to  eliminate  the  excessive  duplication  of  trajectory  programs  and,  more 
importantly,  to  eliminate  the  excessive  time  consumed  by  the  users  in  the 
maintenance  of  familiarity  with  each  of  the  different  programs. 

This  modular  program  consists  of  two  main  sections.  The  first,  containing 
the  executive  routines,  provides  all  of  the  control  logic  for  the  program  from 
the  specification  of  input  data  clear  through  to  the  final  calculated  results. 
Standardized,  general  formats  are  provided  for  the  inclusion  of  all  data.  All 
necessary  standard  mathematical  operations  are  coded  and  Included,  Including 
means  for  the  numerical  integration  of  up  to  28  differential  equations,  as  well 
as  standard  generalized  formats  for  the  printing  of  the  trajectory  results. 

The  second  section  is  written  so  that  each  user  can  select  or  program 
individual  modules  that  meet  his  particular  vehicle  requirements.  The  program 
has  been  written  in  such  a manner  that  it  can  be  used  to  simulate  any  type 
of  flight  in  the  atmosphere  — including  the  simulation  of  guided  vehicles  from 
simple  3 DOF  particle  trajectories  to  maneuvering  6 DOF  simulations  of  air-to- 
air  missiles  with  proportional  navigation  or  maneuvering  re-entry  bodies  flying 
along  evasive  trajectories.  Several  basic  modules  as  well  as  some  specific 
modules  have  been  written  and  are  included  in  this  report  for  the  aid  of  the 
user.  This  program  is  efficient  and  easily  modified  by  the  user  so  that  he 
can  use  it  from  the  original  conception  of  a system  and  its  preliminary  design 
through  to  its  final  flight  evaluation. 

EXECUTIVE  ROUTINES 

The  primary  functions  of  the  executive  routines  are  to  control  the  program 
flow,  establish  standardized  formats  for  the  insertion  of  data  and  modules  into 
the  program,  and  provide  for  the  economic  storage  of  parameters  and  their  use. 

The  flow  logic  of  the  program  is  shown  in  Figure  1.  In  order  to  minimize 
storage,  the  program  was  coded  using  overlays.  The  individual  subroutines 
included  in  each  overlay  level  are  shown  in  Figure  2 and  described  in  the  follow- 
ing section. 

DESCRIPTION  OF  ROUTINES. 

Program  OV.  This  is  the  main,  zero  order  overlay.  Its  function  is  to  call 
the  primary  overlays;  therefore,  it  controls  the  main  flow  of  the  program. 

Subroutine  ZERO.  This  routine  zeros  out  the  entire  common  storage  array, 

Y(l)  through  Y(4940),  and  sets  the  following  default  values: 
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OVERLAY  (TRAJ,  0.  0) 

PROGRAM  OV 

SUBROUTINE  ZERO 

SUBROUTINE  INPUT 

SUBROUTINE  RESET 

SUBROUTINE  TAPTAB 

BLANK  COMMON 

OVERLAY  (TRAJ,  1.0) 

OVERLAY  (TRAJ,  2, 0) 

PROGRAM  START 

PROGRAM  TRJPLTS 

SUBROUTINE  SENCOS 

SUBROUTINE  MYPLOT 

SUBROUTINE  FNOL3 

SUBROUTINE  ARKTAN 

SUBROUTINE  MATINV 

SUBROUTINE  ARDCFT 

SUBROUTINE  FRMRAN 

SUBROUTINE  MATVEC 

SUBROUTINE  ITAB 

OVERLAY  (TRAJ,  1.11 

OVERLAY  (TRAJ.  1.2) 

PROGRAM  SETUP 

PROGRAM  INTGRT 

SUBROUTINE  IMOD1 

SUBROUTINE  DERIV 

• 

SUBROUTINE  TERM 

• 

SUBROUTINE  OUT 

• 

SUBROUTINE  MODI 

SUBROUTINE  IMOD20 

• 

• 

SUBROUTINE  MOD20 

SUBROUTINE  PROCESS 

FIGURE  2 OVERLAY  FILES 
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Y(2302)  = J = 2 
Y(2304)  = XNE  = 0 
Y(2305)  = MPR  = 1 
Y(2306)  = ERROR  = 1 


Integration 

Controls 


Y(3000)  = Rg  = 20925631.  ft  = earth's  radius 

Y(3001)  = W = 7.29211508x10  ^ rad/sec  “ earth's  rotational  velocity 

ti 


These  will  be  the  values  for  the  array  Y(l)  through  Y(4940)  unless  they  are  set 
differently  in  the  data. 


Subroutine  INPUT.  All  data  cards  for  the  initial  run  or  stage  are  read  and 
Interpreted  in  this  subroutine  (see  section  on  Control  Cards). 


Subroutine  RESET.  All  additional  data  for  stages  or  runs  other  than  the 
initial  one  are  read  here. 


Subroutine  TAPTAB.  This  routine  reads  the  tabulated  data  from  UNIT  5 (cards), 
arranges  them  into  an  array,  and  then  BUFFERS  them  out  onto  UNIT  9.  (A  description 
of  how  to  set  up  the  tables  is  included  in  the  section  on  Control  Cards,  page  13. ) 

Program  START.  This  routine  sets  all  of  the  initial  conditions  for  starting 
the  numerical  integration. 

Subroutine  SENCOS  (A,  SA,  CA,  N) . This  subroutine  supplies  sin  A = SA  and 
cos  A = CA  for  the  range  0 ^ A _<  360.  If  N = 0,  A must  be  in  degrees  and  if 
N = 1,  A must  be  in  radians.  Subroutines  SIN  and  COS  must  be  included  as  part 
of  the  system  library. 

Subroutine  FN0L3  (J,  NN,  G,  L,  MPR,  XNE,  T,  C,  D,  DERIV,  TERM,  OUT) . This 
subroutine  (Reference  1)  numerically  integrates  all  of  the  differential  equations. 
The  items  which  need  to  be  defined  are  in  the  following  storage  locations: 

J=Y(2302) 

G=Y(2300) 

MPR  = Y(2305) 

ERROR=Y(2306) 

XNE  = Y(2301) 

These  parameters  are  defined  as: 

J:  (INPUT,  INTEGER) 

This  parameter  indicates  the  integration  method. 


^Ferguson,  R.  E.  and  Orlow,  T.  A.,  "FN0L3,  A Computer  Program  to  Solve  Ordinary 
Differential  Equations,"  NOLTR  71-2,  1 Mar  1971 
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Use  Runge-Kutta  method  of  integration  to  termination. 

Truncation  errors  are  not  calculated;  the  step  size  G is 
not  adjustable. 

Use  Runge-Kutta  for  the  first  three  steps,  then  Adams-Moulton 
for  the  remainder  of  the  interval  of  Integration.  Truncation 
errors  are  calculated.  The  step  size  is  adjustable  unless 
XNE  = 0.  If  the  step  size  is  adjusted,  new  starting  values 
are  obtained  through  the  Runge-Kutta  method. 

Use  Runge-Kutta  throughout.  The  truncation  errors  are 
calculated;  the  step  size  is  adjustable  unless  XNE  >■  0. 

(INPUT,  REAL) 

This  is  the  initial  step  size. 

(INPUT,  INTEGER) 

This  is  the  print  frequency  — the  number  of  integration 
cycles  between  printouts.  If  MPR  = 0,  then  printing  is 
determined  by  values  assigned  to  Y(2998)  and  Y(2997),  where 
Y(2998)  is  set  equal  to  some  running  variable  like  T,  C(l), 

D(l),  etc. and  Y(2997)  is  a constant  Interval  in  Y(2998) 
between  printing  cycles. 

(INPUT,  REAL) 

This  is  the  step  size  control.  The  step  size  is  unchanged  if 
the  worst  of  all  the  errors  lies  within  the  window  10~XNE-3^ 
]^0“XNE  ^ fhe  step  size  is  increased  if  the  errors  are  all 

less  than  10“^^^”^.  The  step  size  is  decreased  if  for  some 
differential  equation  the  error  is  greater  than  10"^^®. 

If  ERROR  < 0 and  XNE  0. , the  automatic  adjustment  of  the 
step  size  is  a function  of  the  absolute  errors. 

If  ERROR  = 0.  and  XNE  ^ 0.,  the  automatic  adjustment  of  the 
step  size  is  a function  of  the  relative  errors. 

If  ERROR  = e > 0.  and  XNE  ^ 0. , the  automatic  adjustment  of 
the  step  size  is  a function  of  the  relative  errors  where  the 
relative  errors  are  equal  to  the  absolute  errors  divided  by 
the  maximum  ERROR,  |C(I)().  This  option  removes  the  possibility 
of  using  "small"  functional  values  to  compute  relative  error, 
otherwise  this  option  is  identical  to  the  previous  option  and 
is  to  be  preferred  over  it.  If  XNE  = 0.,  the  step  size  G is 
not  adjustable.  The  other  parameters  are  either  set  internally 
in  the  program  or  are  defined  in  the  section  on  Control 
Cards. 
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Subroutine  ARKT^  (A,  B,  C,  N).  This  subroutine  calculates  arctangents 
defined  as  C = tan“^(A/B) . If  N » 0,  C is  in  degrees  and  if  N » 1,  C is  In 
radians.  The  range  of  C is:  -180  £ C £ +180.  If  A “ 0,  B < 0.  C - -180;  and 
if  B > 0.  then  C = +180.  Subroutine  ATAN  must  be  In  the  system  library. 

Subroutine  MATINV  (A,  B,  C) . This  subroutine  computes  the  transpose  (B) 
and  inverse  (C)  of  the  (3,  3)  matrix  (A).  If  the  determinant  of  A is  zero, 
neither  B or  C is  calculated;  instead,  a comment  is  printed  and  control  is 
returned  to  the  calling  program. 

Subroutine  ARDCFT  (H,  P,  T,  D,  C,  G) . The  earth's  atmospheric  properties 
(Reference  2)  are  supplied  by  this  subroutine  up  to  an  altitude  of  10^  feet. 
Entering  with  the  altitude  (H,  ft.)  the  pressure  (P) , temperature  (T) , 
density  (D) , speed  of  sound  (C)  and  acceleration  due  to  gravity  (G)  are  given 
ratioed  to  their  corresponding  sea  level  values. 

Subroutine  FRMRAN  (TABLE,  NUM,  MFNC,  U,  A).  This  is  a linear  Interpolation 
routine  which  extracts  tabulated  data  from  TABLE,  and  then  with  the  NUM  Independent 
variables  U]^,  U2.  ...  it  linearly  interpolates  or  extrapolates  2^^^^  _ 1 

times  and  supplies  the  MFNC  values  of  the  functions  A. 

Subroutine  MATVEC  (A,  B,  C,  N).  Products  of  matrices,  whose  orders  are  (3,  3) 
and  (3,  1)  are  computed  by  this  subroutine.  When; 

N = 0,  C = AB* 

N = 1,  C = A^b* 

N = 2,  C = AB 

N = 3,  C = A^b 

N = 4,  C = ABT 

N = 5,  C = aV 

The  A,  B,  and  C arrays  are  stored,  column-wise,  starting  at  the  left.  The  symbol  * 
indicates  that  B,  in  these  cases,  is  a (3,  1)  array.  In  all  other  cases  A,  B,  C 
are  (3,  3)  arrays. 

Subroutine  ITAB  (NTAB,  N,  U,  V) . This  routine  selects  the  table  designated 
by  NTAB,  which  is  the  numerical  location  of  the  table  in  the  KTAB  array,  and  for 
N independent  variables  of  U calls  FRMRAN  for  the  linear  interpolation  of  the 
function  V. 

Program  SETUP.  This  program  calls  the  initial  modules  IMODl  through  IMOD20 
as  designated  by  the  code  1 control  cards. 


^U.S.  Standard  Atmosphere,  1969  (NASA,  Dec  1969,  Washington,  DC) 
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Subroutines  IMODl  through  IM0D20.  These  dummy  subroutines  are  Included  so 
that  users  can  substitute  their  own  subroutines  for  calculating  any  Initial 
conditions  that  may  be  needed  before  the  numerical  Integration  Is  started. 

Program  INTGRT.  This  program  sets  up  the  Integration  controls  and  calls 
FN0L3  for  the  numerical  integration  of  the  equations  of  motion. 

Subroutine  DERIV.  This  subroutine  calls  the  appropriate  modules,  MODI 
through  MOD20,  as  designated  by  the  code  1 control  cards. 

Subroutine  TERM.  The  termination  conditions  (code  4 control  cards)  are 
checked  and  if  any  one  of  them  has  been  met,  the  Integration  is  terminated. 

Subroutine  OUT.  The  output  is  prepared  and  printed  here  based  on  the 
Information  included  on  the  code  2 control  cards. 

Subroutines  MODI  through  MOD20.  These  are  dummy  subroutines.  The  user 
should  substitute  his  own  subroutines  for  the  dummy  ones.  These  subroutines 
are  to  contain  all  of  the  definitions  for  all  of  the  differential  equations. 

Subroutine  PROCESS.  Again,  this  is  a dummy  subroutine  which  can  be  replaced 
by  the  user.  Any  accessary  calculations  that  are  not  needed  for  the  Integration 
of  the  differential  equations  are  usually  Included  in  this  subroutine.  The 
subroutine  is  called  In  subroutine  OUT  everytlme  that  the  print  conditions 
have  been  met. 

Program  TRJPLTS.  This  program  contains  the  calls  for  plotting  any  of  the 
data  designated  on  the  code  5 control  cards. 

Subroutine  MYPLOT.  This  is  a dummy  subroutine.  If  the  user  wishes  to  plot 
any  variables  he  must  substitute  his  own  MYPLOT  subroutine  containing  his  own 
GOULD  or  equivalent  plot  calls. 

A FORTRAN  listing  of  these  executive  routines  has  been  Included  in  Appendix  A. 

STORAGE  ALLOCATION.  As  mentioned  earlier,  the  program  has  been  coded  using 
overlay  files  in  order  to  minimize  the  machine  memory  required.  The  amount 
of  storage  needed  for  execution  will  vary  according  to  the  size  of  the  particular 
modules  used  as  well  as  by  the  size  of  the  arrays  that  are  to  be  plotted  and 
dimensioned  in  subroutine  MYPLOT.  Generally,  on  WOL's  GDC  6500,  the  user  has 
needed  around  45000^0^  locations.  The  maximum  has  been  on  the  order  of  63000^^^ 

locations. 

All  of  the  parameters  stored  in  the  program  have  been  placed  into  an  array 
dimensioned  Y(4940).  The  first  2299  locations  have  been  allocated  to  trajectory 
parameters.  These  are  available  for  coding  in  the  modules.  The  locations  Y(2300) 
through  Y(4940)  are  utilized  in  the  executive  routines  and  as  such,  are  not 
available  to  the  user.  As  an  aid  to  keeping  track  of  the  parameters,  the  Y array 
has  been  broken  into  several  parts.  It  is  not  absolutely  necessary  for  the  user 
to  retain  this  designation  in  his  modules;  but,  it  is  a great  aid  to  keeping 
all  modules  interchangeable.  This  Y array  breakdown  is  shown  in  Figure  3, 

The  fixed  storage  assignments  are  listed  in  Appendix  B. 


11 


NSWC/WOL  TR  78-59 


NSWC/WOL  TR  78-59 


CONTROL  CARDS.  All  data  necessary  for  running  this  program  are  entered 
Into  the  program  through  the  use  of  standardized  formats.  Each  piece  of 
Information,  except  for  tables.  Is  entered  on  a separate  control  card.  The 
particular  use  for  each  data  card  Is  determined  by  the  code  punched  In  the 
first  two  columns  of  the  card.  The  codes,  formats,  and  types  of  Information 
are  tabulated  In  Table  1 on  page  14  and  explained  In  more  detail  below. 

Code  0.  Each  data  deck  must  contain  this  card  as  the  first  card  In  the 
deck.  The  users  title  located  In  columns  3 through  72  will  appear  In  the 
heading  at  the  top  of  each  page  of  printout. 

Code  1.  These  are  the  modules  to  be  used  for  this  particular  run  or 
stage.  These  cards  must  be  read  In  the  order  you  wish  the  module  to  be  called. 
It  Is  not  necessary  that  the  mod  numbers  be  In  sequence. 

Code  2.  These  are  the  variables  to  be  listed  In  the  output.  The  data  will 
be  listed  In  columns  with  the  first  column  always  being  time.  The  first  IS 
variables  will  be  listed  on  the  first  page  and  the  next  15  on  a separate  page. 
The  number  of  variables  listed  can  be  any  number  from  1 through  30;  but,  note 
that  If  you  print  the  results  of  no  more  than  15  variables  you  will  save  paper, 
time,  and  money.  The  columns  will  be  printed  In  the  order  that  the  code  2 
control  cards  appear  In  the  deck.  Each  code  2 control  card  Is  to  contain 
the  location  of  the  parameter  In  the  Y array,  the  heading  that  you  wish  to 
appear  at  the  top  of  the  column,  and  the  format  of  the  parameter.  The  maximum 
width  of  each  column  Is  8 spaces  but  you  can  place  the  decimal  depending 
on  what  Is  being  listed.  If  the  format  Is  left  off  the  card,  the  default  Is 
F8.0. 


Code  3.  These  cards  contain  the  Initial  values  of  any  parameters  In  the 
program.  They  may  be  in  any  order  but  the  total  number  for  all  stages  must 
not  exceed  200. 

Code  4.  This  code  identifies  the  termination  conditions  for  the  numerical 
Integrations.  The  program  will  stop  whenever  any  parameter  in  the  Y array 
designated  as  a stop  variable  goes  outside  of  the  lower  or  upper  limits  as 
set  by  this  code.  As  many  as  10  termination  conditions  may  be  set  for  any 
complete  trajectory  run. 

Code  5.  Any  variable  In  the  Y array  that  is  to  be  plotted  (other  than  time) 
must  be  designated  with  a code  5 control  card.  Whenever  any  code  5 control 
card  (up  to  a maximum  of  10)  appears  in  the  deck,  subroutine  MYPLOT  will  be 
called  and  the  users  plot  options  will  be  performed. 

Codes  6 and  7.  A maximum  of  28  differential  equations  can  be  designated 
with  these  codes.  The  code  6 designates  the  dependent  variables  and  the  code  7 
their  derivatives.  The  variable  locations  In  the  "C"  and  "D"  arrays  must 
correspond;  i.e.,  the  first  variable  in  the  "D"  array  must  be  the  derivative 
of  the  first  variable  in  the  "C"  array. 


TABLE  1 CONTROL  CARDS 
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STOP,  END  OF  DATA 
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Code  8.  These  control  cards  are  used  to  Indicate  where  the  tabulated 
functions  appear  In  the  deck  of  tables.  When  the  modules  are  coded,  an  Index 
In  the  KTAB  array  Is  assigned  to  each  tabulated  function.  For  example,  maybe 
the  axial  force  coefficient  was  coded  as  having  Index  number  12  In  the  KTAB 
array.  For  this  particular  run,  the  axial  force  coefficient  table  that  you 
wish  to  use  may  be  the  fourth  table  In  your  deck  of  tables;  therefore,  IN 
equals  4 and  VAR  equals  12. 

Code  9.  A card  with  the  number  9 In  column  2 must  proceed  the  deck  of 
tables.  It  causes  the  program  to  read  the  following  cards  as  tabulated  data. 
The  tables  should  be  arranged  as  follows: 


Table  1 


Title  card 
Control  card 

Listing  of  Independent  variables 
Listing  of  values  of  dependent 
function 


(FORMAT  7A10) 
(FORMAT  1415) 
(FORMAT  6E12.7) 

(FORMAT  6E12.7) 


Table  2 


Table  I 
(1=3,  49) 


Title  card 
Control  card 

Listing  of  Independent  variables 
Listing  of  values  of  dependent  function 

repeat  for  as  many  tables 
as  needed 


BLANK  CARD 
BLANK  CARD 


at  end  of  tables 


In  this  program  the  tabulated  functions  are  functions  of  1,  2,  or  3 variables, 
with  each  function  located  In  a separate  table.  The  tabulation  of  a function 
of  three  variables  would  be  as  follows: 

a.  Control  Card 

L,  N,  M,  nj^,  n2,...nj^  (FORMAT  1415)  where: 


L = 

0 

all  cases 

N = 

3 

number  of  Independent  variables 

M = 

1 

each  table  contains  only  1 function 

til. 

^2  * ^3 

numbers  of  values  of  each  Independent  variable 
for  which  values  of  the  function  are  tabulated. 
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b.  Listing  of  values  of  Independent  variables  for  which  function  Is 

tabulated 

On  the  first  card(s)  the  n.  values  of  the  first  Independent  variable 
are  listed.  On  succeeding  cards  the  02  and  n^  values  of  the  second  and  third 
independent  variables,  respectively,  are  listed. 

The  following  restrictions  apply:  Values  of  two  different  independent 
variables  may  not  appear  on  the  same  card.  At  least  two  values  of  each 
Independent  variable  must  be  listed;  and,  all  values  must  be  distinct  and 
must  be  listed  in  ascending  order.  The  format  for  these  cards  Is  6E12.7. 

c.  Listing  of  values  of  dependent  function 

The  three  independent  variables  are  designated  N2,  and  N^.  The 
number  of  values  of  each  of  these  variables  Is  n^^,  02,  and  n.,  respectively. 

A block  of  data  contains  those  values  of  the  function  for  all  values  of 
listed,  and  for  one  particular  value  of  N.  and  N2.  The  first  block  corresponds 
to  the  first  value  of  N.  and  N2  listed;  the  second  block  corresponds  to  the 
first  value  of  N,  and  the  second  N2,  etc.  These  blocks  are  repeated  until 
a set  of  blocks  tor  the  first  value  of  and  all  values  of  N2  have  been 
presented.  Sets  for  the  remaining  values  of  follow  until  the  table  is 
completed.  As  a check  there  are  n.  sets,  n.  x 02  blocks  and  n^^  x n2  x n^ 
distinct  values  of  the  function.  The  format  for  this  tabulation  is  also^6E12.7. 

Code  10.  This  card  is  to  be  placed  at  the  end  of  the  data  for  that 
particular  run  or  stage  of  the  run.  When  this  card  is  read,  the  program  will 
stop  reading  data  cards  and  start  integrating  the  equations  of  motion.  When 
one  of  the  termination  conditions  has  been  met,  the  program  will  stop  integrating 
and  start  reading  the  next  coded  data  card.  At  this  point  in  time,  the  program 
still  retains  the  initial  conditions  as  read  in  at  the  beginning  of  the  run 
as  well  as  all  of  the  values  as  last  calculated  when  the  termination  conditions 
were  met. 

If  you  wish  to  stack  runs,  i.e.,  start  another  run  with  slightly  different 
initial  conditions,  follow  the  code  10  card  with  a new  code  0 title  card,  and 
then  follow  this  with  the  necessary  changes  that  you  wish  to  make  in  the  original 
code  3 data.  The  program  will  retain  the  original  initial  code  3 conditions 
except  for  those  that  you  change  here.  You  must  include  new  code  1,  4,  and  8 

cards;  i.e.,  tell  the  program  which  modules  to  use,  new  stop  conditions,  and  which 
tables  to  use. 

Code  11.  This  is  a title  card  that  is  an  indication  to  the  program  that 
staging  is  to  take  place.  The  general  procedure  is  that  the  program  will  read 
additional  data  at  this  time.  These  data  will  then  replace  the  values  retained 
by  the  program  when  the  last  portion  of  the  flight  was  terminated.  This  allows 
you  to  restart  the  calculations  where  you  finished.  The  only  code  3 data 
required  are  those  that  you  wish  to  change  at  that  point  in  the  trajectory. 

You  must  include  new  code  1,  4,  and  8 cards,  i.e.,  tell  the  program  which  modules 
to  use,  which  tables  to  use,  and  new  stop  conditions.  Remember  though,  the 
total  number  of  data  variables,  code  3 cards,  must  not  exceed  200  and  the  total 
number  of  tables  must  not  exceed  49  for  all  stacked  runs  or  stages. 
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Code  99.  This  stops  the  program. 

TRAJECTORY  MODULES 

The  primary  reason  for  writing  this  program  was  to  build  a program  which 
could  be  utilized  and  changed  by  a wide  variety  of  users  without  them  having 
to  spend  an  inordinate  amount  of  time  in  learning  and  adopting  the  code.  It 
was  also  envisioned  that  the  program  had  to  be  of  use  to  those  conducting 
preliminary  design  studies  (when  only  the  basic  fundamentals  of  the  vehicle 
are  known  and  the  vehicle  characteristics  are  constantly  changing)  as  well 
as  for  those  analyzing  the  flight  mechanics  of  production  systems.  In  order 
to  refrain  from  writing  a general  purpose  program  that  would  cover  as  much 
detail  as  possible  for  all  users,  but  satisfy  none,  it  was  decided  to  program 
the  problem  so  that  different  modules,  written  for  specific  systems,  could 
be  selected  or  written  and  Inserted  by  each  user. 

In  order  to  keep  each  module  as  universal  in  its  use  as  possible,  it  was 
necessary  to  break  the  system  model  into  several  parts  and  to  minimize  the 
linkage  between  each  part.  In  general,  the  parts  of  a system  can  be  divided 
as  is  shown  in  Figure  4. 

The  logic  behind  these  parts  was  as  follows.  The  only  information  that 
would  be  passed  from  the  target  module  to  the  seeker  module  would  be  the 
target  coordinates.  The  Interface  between  the  seeker  and  autopilot  modules 
would  usually  consist  of  a maximum  of  three  error  signals.  The  Interface 
between  the  autopilot  and  the  force  and  moment  modules  would  be  the  two  or 
three  control  deflections  and  the  only  thing  that  would  have  to  be  passed 
over  to  the  equations  of  motion  are  the  three  forces  and  three  moments. 

Note  though  that  this  arbitrary  division  of  the  modules  while  appearing 
to  be  logical  is  not  permanently  locked  into  the  code.  For  any  simulation 
as  many  as  20  modules  can  be  used  and  the  quantities  exchanged  between  modules 
can  be  chosen  at  the  whim  of  the  programmer.  These  are  just  suggested  modules 
which  will  aid  in  the  exchange  of  modules  among  all  of  the  users. 

The  direction  cosine  matrix  is  generally  defined  in  the  section  called 
"Matrix  Transformations."  When  using  this  program  for  6D0F  simulations  the 
position  of  the  body  principal  axes  with  respect  to  the  inertial  axes  is 
generally  defined  by  integrating  the  elements  of  the  direction  cosine  matrix. 
These  elements  and  their  derivatives  are  defined  in  this  module  as  well  as 
other  transformation  matrices  such  as  the  Inertial  to  local  axes  transformation, 
and  the  local  to  principal  axes  transformation.  Generally,  one  of  two  modules 
will  take  care  of  the  matrix  transformations.  One  of  these  is  for  3D0F  and  the 
other  for  6D0F.  These  have  been  coded  and  are  Included  in  this  report  as  MODI 
and  M0D4 . 
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The  target  module  Is  also  of  a general  nature.  It  allows  for  a target 
to  either  remain  fixed  with  respect  to  the  earth  or  to  fly  at  constant  velocity 
over  the  earth.  This  general  module  is  also  included  in  this  report  as  MODS. 

If  the  user  desires  any  particular  form  of  maneuvering  target  he  can  easily 
insert  his  own  module. 

The  seeker  and  autopilot  are  more  system  oriented;  therefore,  they  are 
generally  coded  for  a particular  system.  Some  examples  are  included  in  this 
report  as  M0D6  and  M0D7. 

Likewise,  the  forces  and  moments  generated  on  the  vehicle  are  dependent  on 
what  kind  of  control  system,  etc.,  the  vehicle  has.  It  could  have  canard  controls, 
aft  fins,  swivel  nozzle,  etc.  A sample  has  been  Included  as  MODS. 

The  equations  of  motion  are  of  a more  general  nature;  therefore,  examples 
for  modules  for  both  3D0F  and  6D0F  are  included  in  this  report  as  MODs  14  and  9. 

Remember,  you  only  need  to  call  the  modules  which  pertain  to  your  case.  In 
early  studies  of  a preliminary  design,  you  may  only  need  one  or  two  very  simple 
modules.  As  the  system' is  developed  you  can  then  expand  and  add  to  your  module 
package.  You  can  also  take  advantage  of  another's  modules  that  have  been  developed 
if  all  parties  maintained  the  general  interchangeability  features  as  outlined  in 
this  report. 

AXIS  SYSTEMS  AND  TRANSFORMATIONS.  In  general,  only  four  different  axis 
systems  are  needed.  These  are  inertial,  local,  principal,  and  geometric  axis. 

These  can  be  used  to  specify  orientations  of  the  missile  with  respect  to  the 
target  and  with  respect  to  the  earth. 

Inertial.  The  inertial  axes  are  defined  as  follows.  The  origin  is  at  the 
center  of  the  earth  and  the  Xj^  axis  goes  through  the  north  pole  and  the  Yj^  and 
Zr  axes  go  through  the  equator  as  is  shown  in  Figure  5.  These  axes  are  the 
ones  in  which  the  force  equations  are  Integrated. 

Local.  The  local  axes,  labeled  by  the  subscript  "L"  have  their  origin  on 
the  earth's  surface  right  below  the  vehicle.  The  Zj^  axis  is  perpendicular  to 
a tangent  plane  located  in  the  earth's  surface  and  passes  up  through  the 
vehicle's  center  of  gravity.  The  Xl  axis  lies  in  the  tangent  plane  and  always 
points  north.  The  origin  of  the  local  axes  are  defined  by  the  longitude,  tjj, 
and  latitude,  of  the  vehicle.  These  relationships  are  shown  in  Figure  5. 


The  transformation  matrix  for  transferring  a vector  in  the  inertial  system 
to  one  in  the  local  axes  is  designated  [JirJ  where. 
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['«J 


sinii/j^sinTj^ 


-cosi|/_slnT_ 
K K 


-slnt/<j^C08Tj^ 


COSllij^COSTjj 


The  winds  are  generally  defined  with  respect  to  these  local  axes  as  shown 
In  Figure  8. 


Principal  Axes.  These  are  the  principal  axes  of  the  vehicle,  i.e.,  their 

origin  is  at  the  center  of  mass  and  they  consist  of  a set  of  cartesian  axes 

for  which  the  inertia  tensor  will  be  a diagonal  (see  Figures  6 and  7).  The 

initial  orientation  of  the  principal  axes  with  respect  to  the  local  axes  is 

defined  by  the  three  angles  y,,,  and  (in  that  order)  where, 

MM  M 


■ [‘lp] 


cose^cosv.. 
M M 


-coscj^sinyj^ 


[\p]=  cosi-Msinyj^  - sini>„sinc^cosY^  cos^j^cos^  +slnCj^sinYj^  sln(l>^cosej^ 


-siny^j^sinYj^-  cosv^^  sinCj^  cosYj^ 


-sin4!,,cosYw  + costv-slnc  sinY„  cos4'„cose,, 

MM  M n n M M“ 


The  means  of  arriving  at  this  matrix  after  the  initial  time  is  defined  and 
explained  in  M0D4. 

Geometric  Axes.  These  are  a set  of  orthogonal  axes  which  are  defined  for 
maximum  convenience  in  expressing  the  vehicle  aerodynamics.  They  will  generally 
define  a plane  of  symmetry  for  the  vehicle  external  configuration  and  their 
origin  will  be  located  at  the  moment  reference  center.  The  relationship  between 
these  axes  and  the  principal  axes  are  shown  in  Figure  7.  The  center-of-gravity 
(c.g.)  of  the  vehicle  is  defined  with  respect  to  the  origin  of  the  geometric  axes 
by  the  lateral  transformations  X , Y and  Z . The  angular  orientation  is 
expressed  through  the  rotations  and  These  transformations  are 

performed  in  the  force  and  moment  module,  MODS.  The  transformation  matrix  is 
defined  as. 
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COSljJpCOS  p 
cos4>gSin0g 
-sini{)gSln0g 


sln>f„siniti„cos0 
cos^  slniji  COS0 

(j  Vj  I 


-cos;'^sln  ^ 

cosi}'^,cos0g  + 8ln4i^sinipgSln0g 
-sinij-pCosOp  + cos4igSinil»pSlnej, 


sln^^ 

sinijigcosi/zg 

COSi}i^COS(J/g 


Aerodynamic  coefficients  and  orientations  of  the  vehicle  with  respect  to  the 
flow  are  generally  defined  with  respect  to  these  geometric  axes  as  is  shown  in 
Figure  9. 

THREE-DECREES-OF-FREEDON  MODULES . The  logical  procedure  in  normal  point- 
mass  trajectory  calculations  is  to  integrate  the  body  accelerations  in  an  inertial 
system.  In  the  sample  case  presented  in  this  report,  three  modules  are  used 
to  define  the  second  order  equations  of  motion.  These  are  MODs  1,  2,  and  lA 
or  MODs  1,  3,  and  14.  As  mentioned  earlier,  it  is  sometimes  necessary  to  perform 
some  calculations  prior  to  starting  the  actual  trajectory  calculations.  These 
types  of  calculations  are  usually  necessary  in  order  to  provide  the  initial 
conditions  for  the  differential  equations.  Since  these  need  to  be  made  only 
once,  they  are  programmed  as  IMODs  and  are  called  by  the  executive  routine  prior 
to  actually  starting  the  trajectory  calculations.  For  all  MODs  called  by  the 
program  (identified  by  code  1 cards),  the  program  checks  to  see  if  there  is  a 
corresponding  IMOD.  For  the  3D0F  calculations  (as  presented  in  examples  included 
in  this  report)  an  IMOD  is  required  for  MODs  1,  2,  and  3. 

The  flow  logic  used  in  these  examples  is  shown  in  schematic  form  on  the 
next  page.  Remember,  while  the  schematic  is  general  in  nature,  the  MODs  used 
in  this  report  are  only  typical.  Each  user  can  add  to  them  or  replace  them 
with  those  of  his  own  choosing. 

SIX-DEGREES-OF-FREEDOM  MODULES.  For  6D0F  simulations  there  are  likely  to  be 
six  modules  or  more.  Tf  the  complexity  of  a 6D0F  simulation  is  warranted,  it  is 
usually  necessary  to  prepare  a set  of  modules  which  represent  a specific  missile 
system.  In  general  the  simulation  would  be  laid  out  as  shown  in  Figure  4.  A 
group  of  typical  modules  have  been  selected  and  enclosed  with  this  report  in  order 
to  aid  the  user  in  the  compilations  of  modules  for  his  own  simulation.  The  logic 
of  these  sample  modules  is  shown  on  the  following  page. 

IM0D4  and  M0D4  are  used  to  calculate  the  direction  cosine  matrix.  It  is 
unlikely  that  any  changes  or  additions  would  have  to  be  made  to  these  modules 
since  they  deal  with  the  mechanics  of  the  situation,  not  a particular  piece  of 
hardware.  M0D5  is  used  to  calculate  the  location  of  the  target.  Here  again, 
unless  a particular  maneuver  of  the  target  is  to  be  programmed,  it  is  unlikely 
that  changes  would  be  made  to  this  module.  IM0D9  and  M0D9  are  used  to  calculate 
the  equations  of  motion  of  a vehicle  for  which  the  cross  products  of  Inertia 
are  zero.  M0D19  should  be  used  in  place  of  M0D9  if  the  cross  products  of  inertia 
are  not  zero. 
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PROCESS.  This  routine  is  used  to  calculate  any  quantities  that  are  desired 
for  output  but  which  are  not  necessary  for  the  running  of  the  program.  An  example 
of  the  type  of  things  which  can  be  calculated  here  are  the  longitudinal  range 


along  the  equator,  R , and  the  latitudinal  range,  R 


These  are  calculated  as: 


E 


•"e  ^R 


A calculation  of  the  total  distance  traveled  over  the  earth's  curved  surface 
(projection  of  vehicle's  path  on  the  earth's  surface)  is  calculated  in  the 
following  manner. 


^E  “ ■'r  " V 


» (Rg  sin  - Rj.  sin 


2 .2 

AY  = [-Rg  cos  sin  + R^  cos  sin  Tj.] 

2 2 
AZ  = [Rg  cos  cos  Tg  - Rg  cos  4)^  cos  tg] 


where  4;  and  t are  the  last  latitude  and  longitude  calculated  (last  time  program 
executed  subroutine  PROCESS).  The  increment  of  the  chord  of  travel  over  a segment 
of  the  earth's  surface  is  then  expressed  as, 

AC  * /aX^TaY^  + A? 


and  the  angle  subtended  by  the  chord  as 


6 = 2sin 

c 


The  increment  of  surface  traversed  is  then  equal  to 

ARg  = Rge 


and 


^S  = ‘^S  ^*^5 


is  the  total  distance  traveled. 

Any  similar  calculations  can  be  made  and  added  to  PROCESS.  See  Appendix  P 
for  its  current  form. 
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SETUP  PROCEDURE 

The  required  parameters  In  the  setup  deck  will  depend  on  what  modules  are 
In  use;  but,  the  general  procedure  and  types  of  cards  are  usually  the  same. 

The  following  Is  a generalized  listing  of  the  most  used  setup  and  control  cards. 
The  general  procedure  Is  to  prepare  the  deck  as  follows: 

1.  Title  card,  code  0. 

2.  Select  the  appropriate  modules  for  your  simulation  and  arrange 
the  code  1 control  cards  In  the  order  that  the  modules  are  to 
be  called. 

3.  Provide  the  data  cards,  code  3,  required  by  the  executive 
routines.  These  are  mostly  Integration  controls  and  general 
physical  descriptors.  These  are  located  In  the  executive  and 
generalized  Input  sections  of  the  Y array,  locations  Y(2300)  - 
Y(3099).  See  Appendix  B for  specific,  required  parameters. 

4.  Determine  what  the  termination  conditions  are  and  set  their 
values  with  code  4 control  cards. 


5.  Specify  all  of  the  remaining  Initial  conditions  required  by 
the  modules  with  code  3 cards.  These  are  parameters  which 
generally  have  to  do  with  defining  the  Initial  attitude  of 
the  vehicle. 


a.  Locate  the  local  axes  by  giving  the  longitude,  (deg) 
in  Y(3014)  and  the  latitude,  (deg)  in  Y(3015). 

b.  Define  the  initial  altitude,  h(ft)  in  Y(3013). 


c. 


d. 


Locate  the  principal  axes  with  respect  to  the  local  axes 
by  specifying, 

Yj^(deg)  in  Y(2066) 

e„(deg)  in  Y(2067)  IM0D4 

4>M(deg)  in  Y(2068) 

Define  the  Initial  velocity  with  respect  to  the  local 
axes. 


Vj.(fps)  In  Y(610) 
Yg(deg)  in  Y(2209) 
Yjj(deg)  In  Y(2208) 


IM0D2 
MODS  2,  3 


30 


NSWC/WOL  TR  78-59 


e.  Define  the  angular  orientation  of  the  principal  axes  with 
respect  to  the  geometric  axes  that  you  are  using. 

X(.p(ft)  in  Y(3006)  . 

Ycc(ft)  in  Y(3007)  1 

Z(,f;(ft)  in  Y(3008)  > M0D8 

0c(deg)  in  Y(3019)  \ 

(Jjg(deg)  in  Y(3020)  ' 

^-gCdeg)  in  Y(3021) 

f.  Define  all  other  constants  required  by  the  modules  you  are 
using. 

6.  Decide  what  parameters  you  want  listed  in  the  output  and  identify 

them,  their  titles,  and  formats  on  code  2 control  cards. 

7.  List  all  of  the  dependent  variables  on  code  6 control  cards. 

8.  List  all  of  the  derivatives  for  the  above  dependent  variables 

on  code  7 control  cards. 

9.  Supply  the  tables  that  you  are  using.  Define  the  location  (in  your 
deck)  of  each  of  the  table  array  numbers  on  code  8 control  cards 
and  then  place  all  of  the  tables  in  that  specified  order  behind 

a single  code  9 control  card. 
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•OFCk  Tf.1 

C 01/n/7b  1?.S6.«7  JOHN  HOLMES 

OVERLAY (TRaj, 0.0) 

PROGRAM  OV  (INPUT. output. TAPES. TAPEl<JtTAPEV.TAPf6=r>llTPOT) 
C 8/P/77  JOHN  HOI MES 

COMMON  7(4040) 

COMMON/TAB/7(SO) 

EOUI valence ( Y (??1 1 ) .NOPLOT ) 

100  CALL  7ERO 
poo  CALL  INPUT 

300  CALL  overlay (4HTRAJ. 1 .O.bHPECALL) 

IE  (NOPLOT. (;t.o)  go  to  poo 
500  CALL  RESET 
GO  TO  300 

POO  CALL  0VERLaY(4HTRAJ.P.n.6HPFCA(.L> 

GO  TO  500 
1000  STOP 
END 


SUBROUTINE  ZERO 

C 7/P/74  vIOHN  E.  holmes 

COMMON  Y(4940) 

EOUI valence  ( Y (?30?) . J) . ( Y (2305) .Mpw) 
EOUI VALENCE  ( Y (2304) .L ) . (Y (2301) .XNE ) 
E(jUIVALENCF  ( Y (2308)  .EoROP) 

EOUI valence (Y(300n).RE).(Y(300l).WE) 
DO  1 1=1.4040 
1 Y(I)=0.0 

default  options 

J=P 

Y (2302) =P.O 
MPR=l 

Y (2305) =1 .0 
I =0 

Y (2304) =0.0 
XNE=0.0 
ErROR=-1 .0 
WF=0. 000072921 ISOP 

mean  radius  eor  spherical  earth 
RE  = 2092b631  . 

C 

RE  TURN 
END 


SUBROUTINE  INPUT 
C 7/2/74  JOHN  E.  HOLMES 

COMMON  7(4940) 

INTEGER  OUTNO.STPNO  ( 1 0 ) . PI.  OT  ( 1 0 ) . C V AP  ( 3 1 ) . D V AP  ( 3 1 ) 

EOU I valence  (Y  (230  7)  .NOhOD)  . (Y(2312)  .MMOOd)  ) 

EOUI valence ( Y (2300) .NOOUT) . ( Y (2332) .RNMEl ( 1 ) ) . ( Y (23Qp) .OUTE.O(  1 ) ) 
EOUIVALENCE ( Y (2309) . NOT N ) . ( Y ( 24  2? ) . INN0( 1 ) ) . (7(262?) .VAl VE  ( 1 ) ) 
EOU I valence ( Y (231 0) .NOSTOP) .(Y(?0  22).STPNu(l)).(Y(?9  32) .SUP ( 1 ) ) . 
♦ (7(2842)  .Si  0 ( 1) ) 

EOUI valence (Y (2311 ) .NOPLOT ) . (7(2852) .Pi OT ( 1 ) ) 

EOUIVALENCE ( Y (2872) .LOCC ( 1 ) ) . ( Y ( 2903 ) . r V AR ( 1) ) . ( Y (3934) .LUCD( 1 ) ) 
EOUIVALENCE (Y (2965)  .DVAR ( 1)  ) . (7(4890 ) . K T AH ( 1 ) ) . ( Y ( ’87 1 ) .NOTAR) 
EOUIVALENCE (Y(?R7n) .NOnER) . (Y(?H69) .MOVAR) 


TRIOOOIO 

TRlOOlOO 


TR10O140 

TR100150 

TR100160 

TM100170 

TR100180 

TRlbOPOO 

TR100210 

TR100220 

TR100240 
TR100250 
TRIU0260 
TR100270 
TRIC02HO 
TR100290 
TR100300 
TR100310 
TR100320 
TR100330 
TR1003A0 
TR10n350 
TRl003b0 
TRl 00370 
TR10038U 
TR100390 
TR1U0400 
TR100410 
TR100420 
TR100430 
TR1U0440 
TR100450 
TR1U0460 
TR100470 
TR1004R0 
TR100490 
TRIOOSOO 
TRlOOblO 
TR100520 
TR100530 
TR1U0540 
TB100550 
TR100560 
TH100570 
TR1005HO 
TR100590 
TR100600 
TRlUOblO 
TRlU0b20 
TR100630 
TRlOObSO 
TRlUObbO 
TRlOObTO 
TR100680 
TRlOObPO 
TR100700 
TRll)071() 
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tOUl ViLtNCF (V(30S1 ) .K ( 1 ) I . (Y ( 3044) . I I IL  ( i ) 1 
EOUIVALENCF  (Y(a9y4)  «KFNf)) 

EOUIVALENCF (Y (23^?) tHP U 1 ) ) . ( Y ( 2377 ) .Hn2( 1 ) ) 

DIMENSION  TITL (7) .NMOD (20) ,HNMF1 ( 30 ) .OUT  NO ( 30 ) . INNO(200) 

♦ t VALVE  (200)  .Slip  (10)  tSLO(  10)  ,K  (4N)  , 

♦ LOCCOl)  »LOCn(31)  ,KTAH(49) 

DIMENSION  HDl  ( IS)  .HL)2(  )■=>) 

N01N=0 

NOMOOrO 

NOOUT=0 

NnSTOP=0 

NOPLOT=0 

NOVAPsO 

NODER=0 

NOTABsO 

CK=10h 

DO  970  1 = 1. IS 
HDl (I)=5HFO.O 
HD?(I)=5HFO.O 
DO  101  1=1,49 
K(I)=0 

PFAD(S.998)  IP, (TITL ( I ). 1=1 .7) 

X()ATF=DATE  (OHM) 

XTIME=T1ME (HOMY) 

WPTTE (6,5000) 

FORMAT ( IHl ,TS4,*NAVAL  SIIPFACF  WEAPONS  CENTER*/ 

♦ T‘^9,*white  oak  LAHOHATOHY*) 

WRITE (6,5001 ) (T1 TL  ( I ) . 1 = 1 .7) ,XDATE,XTIME 

FORMAT  ( 1HO,10X,7A10,T87.«RUN  [iA1E*,A10,Tl()7.*TIMF*.Alo/) 

FORMATdHl  ) 

FORMAT ( I2,7A1 0) 

IPRT=1 

READ (5, 1000)  IR, TE .VAR.VARR.HEOI ,hEI)?,hE03 
FORMAT! 12, T6.2E1 4.6, A 10, A 7, A 10) 

GO  T0(2,3)  IPRT 
IPRT=? 

1R1  = IR%1N1  = INSVAR1 =VARSVARB1  = VAHR$HFD1 1 =MEn 1 ShFO? 1 =MEft2»HF 03 1 : 

GO  TO  4 

IPRT=1 

WRITE(6,1001) TRl , INI ,VAR1 ,VARP1 ,Ht()) 1 , HE 02 1 , ME  03 1 , I P , I N , V AP , 

♦ VARR,HE01  ,HED2,HEr)3 

FORMAT ( IH  .I?,rA,^E^4.^,A10,A7,Al0,3X,?h*«,3X,I2.IA. 
♦2F14.A,A10.A7,A|0) 

continue 

MOnULF  NUMRFR 

i IF(IR.NE.l)  GO  TO  10 
NOMOD=NOMOn*l 
IF  (NOMOD.6T.20)  GO  TO  6 
GO  TO  7 
WRITE (6,4001 ) 

FORMAT  ( 1H0,*THF  NUMBER  OE  MODULE.S  EXCEEDS  20*) 

STOP 

CONTINUE 
NM0D(N0M0D)=TN 
GO  TO  1 

variable  To  hf  listed  in  output 

1 IF(IR.NE.2)  (?0  TO  15 
N00UT=N00UT,1 
IF (NOOUT.GT.30)  GO  TO  1 1 
GO  TO  12 
WRITE (6,4002) 

FORMAT  ( 1HO.*PPINTOUT  OF  MUWE.  THAN  30  ITEMS  WAS  CONSTDFRED 


Tuiup/ao 
TR100730 
TR10073S 
TR100740 
TH100750 
TR100760 
TR100765 
TH100770 
TH1007B0 
TR1U0790 
TRlUOBOO 
TRIOOHIO 
TR1O0820 
TH]UflB30 
TR100b40 
TR100e42 
TR100844 
TRl 0084S 
TR100U4N 
TR100B50 
TH1O0860 
TR100870 
TR10087? 
TR100874 
TR10087#, 
TR100878 
TR100H79 
TR1008HO 
TRiooas? 
TR100890 
TR1U0900 
TR10090S 
TR100910 
TM10091S 
TH100920 
TH100925 
rMF03TR100930 
TR100940 
TR100945 
TR100950 
TR100955 
TR100960 
TR10096S 
TR10096H 
TR1O0970 
TR100980 
TH100990 
TRIUIOOO 
TRIOIOIO 
TH101()20 
TR101030 
TR101040 
TR101050 
TR101060 
TH101070 
TRIUI 080 
TRIO  1090 
TRIOI 100 
TRIOI 1 10 
TRIO) 120 
TRIOI 130 
TH101140 
TRIOI 150 
TR101160 
TRIOI 170 
TRIOI 180 
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♦EXCESSIVE  4N0  THEY  WERF  OROPPEO*) 

TR10119U 

NOOUT=NOOUT-l 

TR101195 

60  TO  1 

TR1UI200 

12 

continue 

TR101210 

RNNEl (N00UT)=HE01 

TR101220 

IF (HEn3.EQ.CK)  60  TO  230 

TR101222 

IF<N00UT.6T.15)  60  TO  220 

TR101225 

HOI (N00UT)=HF03 

TRI0I227 

60  TO  230 

TR101229 

220 

N0T=N00UT-15 

TR101231 

HD2(NOT)=HFn3 

TR10I233 

230 

CONTINUE 

TR101235 

OUTNO(NOOUT)=IN 

TR101240 

60  TO  1 

TR101250 

C 

TR101260 

C 

DATA  location  AND  VALUE 

TR101270 

C 

TR10I280 

15 

IFtIR.NE.3)  60  TO  20 

TH101290 

NOtN=NOIN*l 

TR1UI300 

IF (NOIN.6T.200)  fiO  TO  16 

TR101310 

60  TO  17 

TR1U1320 

16 

WRITE (6t4003) 

TR101330 

4003 

FORMAT ( 1H0.*N()MRE»  OF  INPUT  VARIABLES  EXCEEDS 

200*) 

TR101340 

STOP 

TR101350 

17 

continue 

TR1U1360 

INN0(N0IN)=IN 

TR101370 

VALVE (NOIN) =VAR 

TRl013e0 

60  TO  1 

TR101390 

C 

TR1U1400 

C 

TERMINATION  VARIABLE  WITH  UPPER  ANO  LOWER  ROUNDS 

TR101410 

C 

TR101420 

20 

IF(IR.NE.4)  GO  TO  30 

TR10143U 

NOSTOP=NOSTOP*l 

TR101440 

IF (NOSTOP. RT. 10)  GO  TO  21 

TR101450 

GO  TO  22 

TR101460 

21 

WPITE (6*4004) 

TR101470 

4004 

FORMAT (IHO. ‘NUMBER  OF  STOP  CONDITIONS  EXCEEDS 

10*) 

TR101480 

STOP 

TR101490 

22 

CONTINUE 

TR101500 

STPNO(NOSTnP)=lN 

TR101510 

SUP (NOSTOP) =VARR 

TR101S20 

St  OtNOSTOP)=VAR 

TR101b30 

GO  TO  1 

TR101540 

C 

TR101550 

C 

variables  to  he  plotted 

TR101560 

c 

TH101570 

30 

IFdR.NE.S)  GO  TO  40 

TR101S80 

NOPLOT=NOPLOT*l 

TR101590 

IF (NOPLOT. GT. 10)  GO  TO  31 

TH101600 

GO  TO  32 

TR1U1610 

31 

WRITE (6*4005) 

TR101620 

4005 

FORMAT (1H0.*NUMPFR  OF  PLOT  VARIABLES  EXCEEDS 

10*) 

TR101630 

STOP 

TR101640 

32 

CONTINUE 

TR101650 

PLOT (NOPLOT) =IN 

TR1U1660 

GO  TO  1 

TR101670 

C 

TR101680 

C 

DEPENDENT  VARIABLES*  C ARRAY 

TH101690 

C 

TRlOl 700 

40 

IF  ( IR.NE.6)  GO  TO  50 

TR101710 

NOVAR=N0VAP*l 

TR101720 

IF (NOVAR.GT.2fl)  GO  TO  42 

TR101730 

GO  TO  43 

TR101740 

42 

WRITE (6*4000) 

TRlOl 750 

4000 

FORMAT ( 1H0.*THE  NUMBER  OF  DIFFERENTIAL  EQUATIONS  txrFEDS  2fl*) 

TH101760 

STOP 

TR101770 

A-4 
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CONTINUE 
LOCC(NOVAB) =VAR 
CVAR(NOVAP>=IN 
GO  TO  1 

C 

C DERIVATIVES,  0 ARRAY 

C 

50  lEdR.NE.T)  GO  TO  60 
NODER=NO0ER*l 
LOCO(NOOER)=VAR 
OVAR(NOOER)=IN 
GO  TO  1 

C 

C TABULATED  VALUES 

60  IF(IR.NE.B)  go  to  70 
NOTAB=NOTAP«I 
IF (NOTAB. GT. 49)  GO  TO  61 
60  TO  62 

61  »<RITE(6,4006) 

6006  FORMAT (lH0,*THe  NUMBER  OF  TABLES  EXCEEDS  4<>*) 
STOP 

62  continue 

rvAR=VAR 
KTAB(IVAR>=1N 
GO  TO  I 

C 

C an  INOICATFR  that  TABLES  ARE  TO  BE  fiFAO 

C 

70  IFdR.NE.S)  GO  TO  80 
CALL  TARTAR  (K  ,KFNI1) 

60  TO  1 

80  CONTINUE 

c end  OF  data  for  a SINGI E Pun 

90  WRITE (6,2000) 

RETURN 

END 


SliPROUTINE  RFSFT 
C 7/2/74  JOHN  E.  HOLMES 

COMMON  Y(4940) 
integer  STPNO(IO) 

EOU I VALENCE (Y (2863) .STOP) 

EOUIVALENCE (Y (2307) .NOmOD) , (Y(P312) ,NM0D(1) ) 

EQUIVALENCE  (Y  ( 2 3 09  ) , NO  TN ) , ( Y ( 2422 ) , INNOd  ) ) , ( Y (2622)  , VAL  VF  (1  ) ) 
equivalence (Y (2871 ) .NOTAB) , (7(4890) .KT ABd ) ) 

EQUIVALENCE (Y(2306) .ERROR) , (Y(2996) .STAGE) 
equivalence (Y (3044) .TITL  d ) ) . ( Y (2301 ) .XNE) 
equivalence (Y(2310) .NOSTOP) 

equivalence (Y (28??) .STPNO( 1 ) ) , ( Y (2832) , SUP ( 1 ) ) . ( Y ( 28  42 ) , SI  0 (1 ) ) 
DIMENSION  TITL  (7)  .NMOr)(20)  .INNl)(200)  , VALVE  (200). KTAD(  49) 
DIMENSION  SUPdO)  .SLOdO) 

READ(5.998)  IP . ( T I TL d ) . I = 1 . 7 ) 

C IP  = 99,  STOP, ALL  PUNS  FINISHED 

IF(IR.EQ.99)  STOP 

C IP  = 11, staging  has  OCCUPPED.  new  moos  and  tables  ape  to  PE  USED 
IFdR.EO.ll)  GO  TO  80 
DO  10  1=1,2307 
10  Y(t)=0.0 

DO  12  1=2862.2868 
1?  Yd)=0.0 
Y(299fl)=0.n 
Y (2999) =0.0 
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TR101780 
TR101790 
TRIOIBOO 
7X101810 
TX10IB20 
TXI01630 
TP101840 
TP10I850 
TR101860 
TP! 01 870 
TR101880 
TR101890 
TR101900 
TR101910 
TR10I920 
TP101930 
TR101940 
TR101950 
TR101960 
TP101970 
TM101960 
TR101990 
TR102000 
TK102010 
TR102020 
TR102030 
TR102040 
TR1020S0 
TR102060 
TH102070 
TR102080 
TP102090 
TP102100 
TR1021  10 
TH102120 
TR102130 
TR102140 
TR102150 
TR102160 
TM102170 
TH102180 
TR102190 
TH102200 
TR102210 
TH102220 
TRl 02230 
TR102240 
TX102250 
TR102260 
TH102270 
TR102280 
TR102290 
TR102300 
TH102310 
TR102320 
TH102330 
TR102340 
TR102350 
TR102360 
TR102370 
TH1023H0 
TR102390 
TR102400 
TH102410 
TR102420 
TP102430 
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REMIND  IV 

“ 

TR102435 

DO  U 1=3002.3043 

TR102440 

u 

Y(I)=0.0 

TRlb245U 

c****< 

TR102460 

c 

DEFAUI  T OPTIONS 

TR102470 

c 

TR102480 

c 

J=2 

TR1UP490 

Y(2302)=2.0 

TR102500 

c 

MPR=1 

TB102510 

Y(?30S»=1.0 

TH102520 

c 

L = 0 

TR102530 

Y(2304)=0.0 

TR1U2540 

XNE=0.0 

TR102550 

ERROR=-l .0 

TR1U2560 

c****< 

TR102570 

GO  TO  85 

TR102580 

80 

continue 

TR102590 

STOP=0. 

TR102b0O 

stage=stagf.i.o 

TR102610 

85 

continue 

TR102620 

WRITE (6.999)  (TITl 

(I) .1=1.7) 

TR102630 

NOMODsO 

TH102640 

NOSTOP=0 

TR1026S0 

NnTA8=0 

TR102660 

100 

RFAD<5.1000)  IR.IN 

.VAR.VARR, 

HEOl ,ME02.hE03. 

HF04 

TR102670 

WRITE (6. lOOl ) IR, IN. VAR.VARR 

,HED1 .HED2. 

HEl)3 

,HFD4 

TR102680 

1000 

F0RMAT(I2.T6.2F14. 

6.Ain.A7,A10.A9) 

TR102690 

1001 

format (IX, T?, 16, ?F 14.6,A10,A7,A10,  A9) 

TR102700 

C 

MODULE  NUMRFR 

TR102710 

IF(IR.NE.l)  GO  TO 

110 

TR102720 

NOMOD=NOMOn*l 

TR102730 

NMOD(NOMOO)=IN 

TH102740 

GO  TO  100 

TR102750 

C 

DATA  LOCATION  AND 

VALUE 

TP102760 

no 

IF(IR.NE.3)  GO  TO 

120 

TR102770 

NOIN=NOIN*1 

TR102780 

INN0(N0IN)=IN 

TR102790 

VALVE (NOIN) =VAR 

TRlOPbOO 

GO  TO  100 

TR102810 

1?0 

IF(IR.NE.4)  go  to 

125 

TR102820 

c 

termination  variables 

TR102830 

NOST0P=N0STOP*l 

TRlb2840 

STPNO (NOSTOP) =IN 

TRiopeso 

SUP (NOSTOP) =VARR 

TR102860 

SL  O(NOSTOP) =VAR 

TR102e70 

GO  TO  400 

TH102880 

c 

TABULATED  values 

TR102890 

125 

IFdR.NE.B)  GO  TO 

130 

TR102900 

N0TAB=N0TAB+1 

TR1U2910 

IVAR=VAR 

TH102920 

KTAB(IVAR)=IN 

TR102930 

GO  TO  100 

TR102940 

c 

END  OF  data  for  a- 

SINCM  E run 

TH102950 

130 

IF(IR.NE.ln)  GO  TO 

100 

TR102960 

WRITE (6. 2000) 

TR102970 

2000 

FORMAT(lHl) 

TR102980 

990 

FORMAT ( 1H0.7A1  0) 

TR102990 

998 

FORMAT(I2.7A10) 

TH103000 

RETURN 

TR103010 

END 

TH103020 

SUBROUTINE  T APTAR ( K . KEND ) 

TR103030 

C 

TR103040 

C 

7/2/74  lOHN  E. 

HOLMES 

TR103050 

C 

TR103060 

C 

IT  READS  The  DATA 

CARDS  for 

THE  tables 

. arranges  THEM  INTO  THE 

TR103070 

C 

TABLE  array.  AND  WRITES  THE 

TABLE  ARRAY 

ON 

Ell  E 9 

TR1030Bn 
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Tartar  is  a suRROUTlNf  to  hi  user  with  krmwan. 


table 

table 

• • • 

table 
table  I 


(n  = N(ll  NtiMHtR  OE  ARGUMENTS  E OB  E?PST  VARIABLE 

(?)  = N(?)  NiJMHfB  OE  AHGUMENTS  E OB  SETONO  VARIABLE 


table  (M)  = N(N)  NUMBER  OE  ARGUMENTS  FOB  N-TH  VABIABLE 

table  <N*n  = U(l*l)  FIRST  OE  ARGUMENTS  C()RRE  SPRNO I NG  TO 

The  EIBST  VABIA8IE 

• • • 

TABLE(N*N(1)  ) = Ud.Nd))  LAST  OF  ARGUMENTS  COPRESPONO ING 

To  the  EIWST  VARIABl.F 

table  (N*Nd  ) *1  1 = U(8d)  FIRST  OE  THE  ARGUMENTS  CORRESPONDING 

TO  The  second  VARIABIE 

• • • 

table (N*Nd ) ♦N(?> ) = U(?,N(2))  1 AST  OF  THE  ARGUMENTS 

CORRESPONDING  TO  THE  SECOND  VAHIABLE 

• • • 

table (N*N d ) ♦NI?) ♦... ♦N (N) » = U(N«N(N)>  I AST  OF  ARGUMENTS 

CORRESPONDING  TO  THE  N TH  VARIABiF 

table  (N»Nd  ) ♦N(?.*.  . .♦N  (N)  ♦ 1 ) = ARGd»l,....I)  TABLE  VALUE 

corresponding  TO  THE  FIRST  ARGUMENT  Of 
EACH  VARIABLE 

TABLE  <N  + Nd  ) ♦N(2>  ♦...  ♦N(N)  ♦!)  = APG  d * J . . . . . 1 » 2 ) 


table  (N*N  d ) ♦N(?)  ♦ , . .♦N  (N)  ♦Nd  ) ♦ I ) = ARG  (1  . I ♦ . . . • 2 . 1) 

TABLE  (N  + N d ) tN  (2)  . .+N(N)  ♦Nd  ) «N(2)  . .*N(N)  ) = ARG(Nd»fN(?) 

• . . . « N ( N ) ) 

READS  In  THE  ADDITIONAL  TAHLE  VALUES  STARTING  IN  THE  NEXT 
available  location  in  TmE  table  array  ANO  using  THE  SAME 
ARRANGEMENT  AS  ABOVE 


FOR  additional  tables*  TARTAb  READS  THEM  IN  AS  ABOVE*  STARTING 
IN  THE  NEXT  available  IOCA1ION  IN  THE  TABLE  ARRAY 


L *NUM*MFNr*N d ) *N (?)*.., *N (NUM)  (WITH  FORMAT  141b)  T 
THE  CONTROL  CARDS  MUST  HE  OF  THF:  EOI  LOWING  FORMAT  T 

7 

M(J)  = number  of  functions  TABULATE!.  FuR  ThE  J-TH  TAHLE  1 
NI(J)=  number  of  VARIABLES  IN  T ML  J-TB  TABI  E T 
L=n*NUM  NF  0 one  TABLE  OE  NUM  InDEPENOENI  VAbIARiES  IS  REAO  ^ 

THE  DECK  OF  THE  TABLE  TS  AS  FOLLOWS  T 
Nd)  ARGUMENTS  (VA(UES)  OF  THf  FIRST  T 
variable  - IN  ASCENDING  ORDEO  WITH  SIX  PER  1 

card  T 
N(?)  arguments  of  the  SFCONf*  VARIABl.t*  1 
beginning  on  a new  card*  in  ascending  OPDFRI 
with  SIX  PER  CARO  1 
• • ■ I 
N(N|IM)  arguments  OF  U(NUM)*  beginning  ON  a I 
NEW  CARO*  IN  ASCENDING  ORDER  WITH  SIX  PEW  I 
CARD.  1 
THE  FIRSI  FUNCTION  TABULATFn  IS  PUNCHED  IN  I 
BLOCKS  OE  n(NUM)  WITH  SIX  PFP  CARD.  ThF  I 
FIRST  block  STARTS  WITH  THF  VALUE  1 
COPRF SPONDING  to  the  SMALLEST  ARGUMENT  OF  ^ 
each  variable  AMD  Allows  U(M1  TO  VARY.  1 
THE  NEXT  BLOCK  STARTS  ON  A NEW  CARD  AND  hASI 
U(N-l)  TnCRFMFNTED  And  allows  U(N)  TO  VARY  ^ 
AND  SO  ON  UNTIL  T BE  ENTIBE  TABULATION  IS  1 
COMPI ETFD.  I 


1U30N0 

103100 

103110 

103120 

103130 

103140 

1U31S0 

lU31bO 

103170 

103180 

103190 

103200 

103210 

103220 

103230 

103240 

1032S0 

103260 

103270 

1U32B0 

103290 

103300 

103310 

103320 

103330 

103340 

103360 

103350 

103370 

103380 

103390 

103400 

103410 

103420 

103430 

103440 

103450 

103460 

103470 

103480 

103530 

103520 

103510 

103500 

103490 

103540 

103550 

103560 

103570 

103580 

103590 

103600 

103610 

103620 

103630 

103640 

103650 

103660 

103670 

103680 

103690 

103700 

103710 

103720 

103730 

103740 


THIS  PASS  15  BSST  QUAIiITI  PRAOTIOiBIil 

FROM  COrY  XLPtlAlSHr.J  TO  DDC 
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IF  MfNC  IS  GRFATfR  THAN  ONE t REPEAT  FOR 
EACH  SUCCEEDING  FUNCTION  TO  BE  INCLUDED 
THE  TABLE. 


TR 
IN  TR 


L=0.NUM=0 


TERMINATE  the  run 


DIMENSION  K(1  ) .TAHLEdnO)  .Ndb)  .M(49)  .T(7)  .NI  (49) 
equivalence  (NARG.AHG) 

KOUNT  = 1 
NFXLOC  = 0 
DO  310  1=1,100 
310  TARLEd)=0. 

300  READ(S,201)  ( T ( 1)  . I = 1 , 7 ) , DllM 

WRITE (6t20?) 

WRITE (G, 201 ) (T (I ) ,1=1 .7) ,DUM 

READ  (S,  100)  L fNUM.M  (KOIiNT  ) . (N  (I  ) , 1 = 1 ,F>IIM)  ,DUM 
WRITE  (6,1  00)  L,NUM,M( KOUNT)  , (Nd)  ,I  = 1,NUM)  ,DUh 
IF  (L.NE.O.OP.NUM.GT.m.OP.M(KOUNT)  .GT.15)  GOTO  800 
IE  (NUM.EO.O)  RETURN 
NT (KOUNT)  = NUM 
K (KOUNT) =NEXI  OC*l 
C 

C TO  CALCULATE  THE  PRODUCT  OF  THE  Nd)  S 

C AND  PUT  THE  Nd)  S IN  THE  TABLE 

C 

N«;UM  = 0 
NPROD  = 1 
DO  400  1=1. NUM 
NARG  = Nd) 

NSUM  = NSUM  ♦ nAPG 
NPROO  = NPPOn  * NARG 
TABLEd)  = APR 
400  CONTINUE 

BUFFER  OUT  ( R , 1 )( TABL E ( 1 ), TABLE (NUM ) ) 

IE(UNIT(9))  450,  450,1000 
450  CONTINUE 


TO  PUT  Ud) 


IN  THE  T ABl E 


KSTART  = NUM  + NFXLOC  + 1 

DO  500  1=1, NUM 

KFND  = KSTAPT  ♦ Nd)-  1 

KKENn=KENO-KSTART,l 

IE  (KKEND.GF.99)  GO  TO  hOO 

RE  AO (5, 1 10)  (TABLE (KK ) ,KK  = 1 .KKEND) 

WPITE (6,201)  (TABLE (KK ) .KK=1 ,KKEND) 

BUFFER  OUT  (9,1) (TABLt (1) ,TABLE(KKEND) ) 

IF  (UNIT(9n  490,490,1000 
CONT INUE 

KSTART  = KFND  ♦ 1 

continue  1 

I 

TO  PUT  THF  TABLUATEO  VALUES  IN  THE  TABLE 

KADD  = N(Num)  - 1 

IFND  = NPRoO  * M(KOUNT)  / N(NUM) 

DO  600  I=l,IENO 

KFND  = KSTART  ♦ KADD 

KKEND=KEND-KSTART,1 

IF (KKEND. GF. OR)  GO  TO  ROO 

READ (5, 1 10)  (TABLE (KK) ,KK=1 , KKEND) 

WRITE (6,203)  ( TABLE  (KK) ,KK  = 1, KKEND) 

BUFFER  OUT  ( 9 , 1 ) ( T ABLE (1 ) , T ABl E ( KKFNO ) ) 
IF(UNIT(9))  580,580,1000 


103750 

103760 

103770 

103780 

103790 

1U3800 

103810 

103820 

103830 

103840 

103850 

103860 

103870 

103680 

103890 

103900 

103910 

103920 

103930 

103940 

103950 

103960 

103970 

103980 

103990 

104000 

104010 

104020 

104030 

104040 

104050 

104060 

104070 

104090 

104100 

104110 

104120 

104130 

104140 

104150 

104160 

104170 

104180 

104190 

104200 

104210 

104220 

104230 

104240 

104250 

104260 

104270 

104280 

104290 

104300 

104310 

104320 

104330 

104340 

104350 

104360 

104370 

104380 

104390 

104400 

104410 


i 
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580 

continue 

TR104420 

KSTART  = KFNO  ♦ 1 

TR104430 

600 

continue 

TR104440 

NFXLOr=KENr> 

TR104450 

KOUNT  = KOmnT  ♦ 1 

TR104460 

IF  (KENO.GE. 3000)00  TO  950 

TH1U4470 

GO  TO  300 

TR104480 

800 

WRITE (6t200) 

TR104490 

STOP 

TR104500 

900 

WRITE (6t220) 

TR1045ln 

STOP 

TR104520 

950 

WRITE (6.230) 

TR104530 

STOP 

TRI04540 

1000 

WRITE (6.240) 

TR104550 

STOP 

TR104560 

100 

FORMAT(l4lc;,?x,ARl 

TH104570 

no 

FORMAT  (6E12.5) 

TR10458n 

200 

FOPMAT(1HO.«  TARTAR  REOUIRES  THAT  L BE  EQUAL  TO  0 AMO  ThAT  N 

TH104590 

♦AND  M BE  BETWEEN  0 AND  15  — ALSO  CHECK  FOR  INCORRECT  NUMBER 

TR104600 

♦OE  DATA  CAPOS  *) 

TR1U4610 

201 

FORMAT  (1X.7A10.1X.A8) 

TR104620 

20? 

FORMAT  (IHO) 

TR1U4630 

220 

FORMATdHO.*  TAPTAB  IS  ATTEMPTING  TO  REAf)  MORE  THAN  V9  VAl  UES 
1AT  ONE  TIMF.  the  DIMENSION  OF  THE  TABLE  ARRAY  WILL  HAVE  TO 

2RF  increased*) 

TR104650 

203 

FORMAT  (lPE.F?n.6) 

TR104640 

230 

FORMATdHO.*  TAHlFS  EXCEED  DIMENSIONED  SI7E  ME  3000  *) 

TR1046H0 

240 

FORMATdHO.*  PARITY  ERROR  OCCURRED  DURING  THE  WPITTmG 

TR104690 

♦ OF  THF  TABi  FS  ONTO  UNIT  9 *) 

END 

TB104700 

OVERLAY(TRaj.I.O) 

PROGRAM  START 

TR1047?n 

C 

7/2/74  JOHN  F.  HUIMFS 

TR104740 

COMMON  Y(4940) 

TR104750 

COMMON/TABI  /table (3000) 

TR10476n 

FOU I VALENCE  ( Y (?309)  .NOIN)  . ( Y (2996)  .STAGE) 

TR1047r0 

EOUIVALENCF (Y (24??) . INNO (I  ) ) . ( Y (262?)  • VALVE  (I  ) ) . ( Y (?9q6) .kFND) 

TR104780 

DIMENSION  IMNO(?nO) .VALVE (200) 

TH104790 

REWIND  9 

TR104800 

DO  lOol  T=l.NOTN 

TR104M10 

IN=INMO ( I ) 

TR10482n 

VAR=VALVE ( T ) 

TR104830 

Y ( IN) =VAR 

TR104840 

1001 

continue 

TR104eS0 

IE (STAGE. GT. 0.0)  GO  Tu  1005 

TR104860 

L(  EN=n 

TRl04e70 

N=1 

TRiOAeao 

M=3000 

TR10489!) 

1003 

BUFFER  IN  (9.1 ) (TABLE (M) .table  (M)  ) 

TB104900 

IF(UNTT(9))  10.12.11 

TR104910 

10 

LFN=LFNGTH (9) 

TR104920 

LI  EN=LEN^LLFN 

TR104930 

IF  (LLEN.EQ.3n0n)  (;0  TO  12 

TR104940 

N=LLEN^1 

TRl0495n 

GO  TO  1003 

TR104960 

11 

WRITE (6»20n0) 

TR104970 

2000 

FORMATdHO.*  PARITY  ERROR  OCCURRED  WHILE  ATTEMPTING  To 

TR104980 

♦ READ  TABLES  FROM  ,|NIT  9 *) 

TR104990 

STOP 

TRlOSOOn 

12 

CONTINUE 

CAL  L OVERLAY (4HTRAJ. 1 . 1 .6HRECALL  ) 

TP105010 

lOOS 

CALI-  OVERLAY(4HTRAJ.1.?.6hrECALI  ) 

END 

TR10G080 

SUBROUTINE  SFNCOS ( A.SA.CA.KEY) 

TRl0509n 

IF(KEY.NE.n)  GO  TO  in 

TH105100 

AA=A/S7.295779S 

A-9 

TRlOSl  10 
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10 


10 

11 

12 

13 

14 

15 

16 
17 
IB 

I 

20 

21 

22 

23 

24 


1 q 


10 


11 


20 

6 


C**** 

C***« 


SA=SIN(AA) 

CA=COS(AA) 

RFTURN 
SA=SIN (A) 

CA=COS(A) 

BFTU9M 

END 

SUBROUTINE  APKTAn  (A,H,C*M0DE) 

IF  ( A ) 1 0 < A . 1 0 
IF  (H)  5.6,6 
7=3.14159265 
fiO  TO  18 
7 = 0. 

GO  TO  21 
IF (B) 13.11.17 
7=SIGN(1. 5707963. A) 

60  TO  18 

Z= A tan (A/B) .SIGN (T. 14169265. A) 

IF (7-3. 1415q?65) 16,15.16 
7 = -Z 

GO  TO  18 
Z=ATAN(A/B) 

IF(MOOE)21.19.21 
0=57.2957795*7 
GO  TO  22 
C = Z 

IF  (ABC  (0-1  .F-0  7 I 23.23.24 
0 = 0.0 
RF  TURN 
FNO 

SUBROUTINE  MATINV(A,R,r) 

DIMENSION  A (9) ,B (9) ,0 (0) 

DFLA=4(l)*4(5)*A(9)*A(2)*A(6)*A(7).4(3)*A(fl)*A(4)- 
1 (P)*A(6)«A(l)-A(9)*A(2t*A(4) 

IF (DFl  A)  10.20.10 
J = 0 

DO  11  1 = 1. 7. 3 
B ( I ) =A ( J*1 ) 

B ( !♦ 1 ) =A ( J.A) 

B(I*2)=A(J.7) 

J=J*1 

0 ( 1 ) = (B (5) *H (9) -H (6) *H (H) ) /DFL A 
0(2)=(B(6)*B(7)-B(4)*H(9) )/DFLA 
0 (3)  = (B (4) *p (8) -B (5) *B ( 7) ) /DEL  A 
0 (4)  = (B (3) *B (B) -B (2) *H (9) ) /DEL  A 
0 (5)  = (B ( 1 ) *P (9) -« (3) *M ( 7)  ) /DEL  A 
0 (6)  = (B (2) *B (7) -B ( 1 ) *H (8) ) /DEL  A 
0 (7)  = (B (2) *H (6) -B (3) *P (5) ) /DEL  A 
0(fl)  = (P (3) *P (4) -B ( 1 ) *H (6) ) /DEL  A 
0(9)=(B(1 )*R(5)-R(?)*B(4) )/PELA 
RFTURN 
WRITE  (6.6) 

F0RMAT(34H  DETERMINANT  OF  A IS  FOUAl  TO  ZERO) 
RFTURN 
FND 

SUBROUT  I NE  A ROOF  T ( H .PPZ  . T TZ  . RRZ  . OCZ  , (.GZ  ) 

DIMENSION  T ( 12) , A ( 1 0) 

ATMOS  OOMPIITFS  STANDARD  PRESSURE.  TEMPERATURE,  AND 
FOR  A GIVEN  altitude  H IN  FEET. 
X(7.A,B.G)=A*4L0G(Z*(Z-H)*0) 

Y(Z.A,P,C)=A*ATAn(7*  P-C> 

DATA ( A ( I ) . 1=1 , 1 0) /-3.35014S769F-1 7.3. 161762924E-14 
1 1 ,2.B48535349F-9.-3.930824139E-7,3.*32295909e-5.-l 
1 j.2564  0363()F-2, -5. 232974b73E-l  ,-3.95524200  7/ 

Z =0.0003048*h 
TM  = 0. 

A-10 


TR105120 

TR105130 

TR10S140 

TR105150 

TR10S160 

TR1U5170 

TR1051RO 

TR105190 

TR105200 

TR1U5210 

TR105220 

TR105230 

TR105240 

TR10525n 

TR105260 

TR105270 

TR105280 

TR105290 

TR105300 

TR105310 

TR105320 

TR105330 

TR105340 

TR105350 

TK105360 

TH10S370 

TR1053B0 

TR105390 

TR10S400 

TR105410 

TR105420 

TR105430 

A (7) *A (5) *A (3) -ATR105440 
TR105450 
TR105460 
TR105470 
TR105480 
TM105490 
TR105500 
TH105510 
TR1U5520 
TR105530 
TR105540 
TRl055bO 
TR105560 
TR105b70 
TR105580 
TR105590 
TR105600 
TR105610 
TR 105620 
TR105630 
TR105640 
TH105650 
TR105660 
TR105670 
TR105680 
DENSITY  RATID5  TR105690 
TR105700 
TR10571O 
TR105720 

. -1. 2699 1 9974F- I TR 105730 

.B3296214bE-3,  TR105740 

TR105750 
TR105760 
TR105770 
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DO  5 1=1*10 
S TM=(TM*A(n )*Z 
Tm=  TM*283.74«)2391 
S1=Z*6356.77 
7(  l)=-l,4655396f-  7/Sl 
T(  2)=  2.5f.S3341E-l  1*AI  06(S1  ) 

T(  3)=  1.4116834E-  4*Al.0G(  14.00238S*Z) 

T(  4)=-3.8282910E-  5*AL0G ( 2 1 6 . 232250-Z ) 

T(  5)=X(2*  1 .50e4978E-4«  26.414270  * 684.10967  ) 

T(  6)=V(Z.  6.74198e0E-4*  0.044294586*  0.5850046  ) 

T(  7)=X(2*  8. 5519675E-5* 137.4745  *10533.544  ) 

T(  8)=y(2*  4.9863416E-S.  0.013120767*  0.90188546J 

T(  9)=X(Z*-?.5392354E-4*193. 32352  *10180.367  ) 

T(10>=Y(Z*  1.1921B79E-7*  0.034567717*  3.3413764  ) 

T(ll)=X(2*-3.3868604E-5*3e4. 32662  *38131.516  ) 

T(12)=Y(Z*  8.9812379E-S.  0.02881021  * 5.5362654  ) 

TINTEG=0.28016067E-02 
DO  10  1=1*12 
10  TINTE6=TINTEG*T< I ) 

(JTNTEG=  -3483.6764*T1nTEG 
TTZ=TM/288. 16 
PPZ=EXP(QImTFG) 

rpz=ppz/ttz 

CrZ=S0PT (401 .874 *TM) /34 0.294 
GG2= ( 1 ./ ( 1 .*2/6356.766) ) **2 
PE  TURN 
END 

SUBROUTINE  FRMRAN  ( T ABLE « NUM . MFNC * U * A ) 

DIMENSION  table ( 1 ) .U ( 1 ) * a ( 1 ) *T ( 1) *TEMP ( 15) *N (15) * IOUM ( 1 ) 
COMMON  /TAR/T 

equivalence  (N ( 1 ) .TEMP ( 1 ) ) , ( IDUM( 1 ) *T ( 1 ) ) 

NUMB  = NUM 

NlNS=2 

NSUM  = 0 

NNPWOO  = 1 

DO  300  1=1. NUMB 

TEMPI  I)  = TABI.Ed) 

NSUM  = NSUM  ♦ N(I) 

NNPROn  = NNPPOO  * N(I) 

300  continue 

to  fill  the  t apray 

NUMT  = 2 * NUMP 
JPOS  = NUMB 
DO  400  1=1. NUMB 

TO  SPACF  TO  THE  BEGINNING  OF  ARGUMENTS  CORPFSPONO I NG  TO 
THE  I-Th  variable 

JSTRT  = JPOS  ♦ 2 
JPOS  = JPOS  ♦ N(I) 

DO  410  J=JSTRT.JP0S 

IF (TABLE ( J) .PT.U ( I ) ) GO  TO  420 

410  continue 

J=JPOS 

4?0  IJ  = NUMT  ♦ I 

inUM(IJ)  = J - JSTRT  *2 

T(l)  THROUGH  T(N)  APE  The  ARGUMENTS  CORRESPONDING  TO  THE  N 
VARIABLES 

7(1)  = TAB(  F(J) 

IJ  = I ♦ NUMB 

A-11 


TR105780 
TR105790 
TR105800 
TR105810 
TH10S820 
TR105830 
TR105840 
TR105850 
TR105860 
TR105870 
TRlOSeBO 
TR105890 
TR105900 
TR105910 
TR105920 
TR105930 
TH10S940 
TR105950 
TR105960 
TR105970 
TR105980 
TR105990 
TR106000 
TR106010 
TR106020 
TR106030 
TR106040 
TR106050 
TR106060 
TR106070 
TR106080 
TR106090 
TR106100 
TR1061 10 
TR106120 
TH106130 
TR106140 
TR106150 
TR106160 
TR106170 
TH106180 
TR106190 
TR106200 
TH106210 
TR106220 
TR106230 
TR106240 
TB106250 
TRl 06260 
TR106270 
TR106280 
TR106290 
TP106300 
TR1063IO 
TR106320 
TRl 06330 
TH106340 
TR106350 
TR106360 
TR106370 
TR106380 
TR106390 
TR106400 
TR106410 
TR106420 
TR106430 
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TtN*l)  through  T(?*M)  ARF  THt  AHGIIMFNTS  ONt  POSITION  HRLOW  THF 
ABOVE  ARGUMENTS 

T<IJ)  = TAHLF(J-n 
AOO  CONTINUE 

ISTART  = NumT  ♦ 1 

NUMTR=NUMT*NUMB 

ISUM=IDUM(NUMTB) -1 

IF  (NUMB.EO.l)  GO  TO  4*0 

NINS=NUMB*1 

NTOP=NUMB-1 

NPROD  = 1 

DO  4T0  I=1.NT0P 

I J=NUMTB-I 

I I J=NINS-I 

NPROn  = NPPOO  * N(IIJ) 

ISUM  = ISUm  ♦ ( I|)IIM  ( IJ) -2)  * NPROO 

430  Continue 

440  TIOUM  = NUMB  ♦ ISUM  ♦ mSUM 
00  1000  M=  l.MFNC 
IDUMI ISTART)  = IinUM 
INDEX  = 1 
NPPOO  = 1 

TO  COMPUTE  the  indites  OF  THF  lAHI.E  VAI  UFS  NEEDED  FOR  THE 
INTERPOl  ATIOM 

00  500  1=1. NUMB 
INDEX  = INDEX  ♦ NIJMT 
DO  510  J=l. INDEX 
IJ  = INDEX  ♦ J 
KJ  = NUMT  ♦ J 

THE  IDUm  array  CONTAINS  THE  VALUES  OF  THE  INOIt'F-*:  OF  THE 
table  values  NEEOEI)  for  THE  INTERPOl  ATION 

lOUM(TJ)  = IDUM(KJ)  ♦ MPHOO 
510  CONTINUE 
I T=NINS-I 

NPROO  = NPROD  * N(II) 

Index  = inofx  •*  index 
500  continue 

TO  PUT  THF  TAHI F VALUES  NEEDED  FOR  THE  INTF RPOL A T ION  IN  THE 
T ARRAY  starting  WITH  T(2*N*1) 

DO  600  I=ISTART.TJ 
KJ  = TOUM(I) 

T(I)  = TABlF(KJ) 

600  CONTINUE 

INTFRPOl  ATION 

JFND  = 2**numB  ♦ ISTART  - ? 

KJ  = NUMB  ♦ 1 
DO  700  1=1. NUMB 
Ij  = KJ  - I 
INDEX  = NUmr  ♦ IJ 

TFM  = (U  ( I J) -T  ( INI'FX)  ) / ( I ( T J) -T  ( INDf  X)  ) 

IJ  = TSTARt 

DO  710  J=ISTART.  JF  Nll.t- 

T(IJ)  = (T ( J. 1 ) -T ( J) ) *TFM  ♦ T(J) 

IJ  = IJ  ♦ 1 


TH106440 
TR106450 
TR106460 
TR106470 
TR106480 
TR106490 
TR106500 
TR106510 
TR106520 
TR106530 
TH106540 
TR106550 
TR1U6560 
TR106570 
TH106580 
TR106590 
TR10b600 
TR1U6610 
TR106620 
TH106630 
TR106640 
TR106650 
TR106660 
TR106670 
TR106680 
TR106690 
TR106700 
TR106710 
TR106720 
TR106730 
TR106740 
TR106750 
TR106760 
TR106770 
TP1067W0 
TR106790 
TRlObSOO 
TR106810 
TRl 06820 
TrI U6830 
TH106840 
TR106850 
TRlObBbO 
TR106870 
TR106880 
TR106890 
TR106900 
TRI06910 
TR106920 
TR106930 
TR106940 
TH106950 
TRl069bO 
TH106970 
TR106980 
TRl 06990 
TRIOTOOO 
TR107010 
TR107020 
TRlOTOJO 
TR107040 
TRl 07050 
TR107060 
TR10707O 
TH1070H0 
TH107090 
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710 

CONTlNUt 

Thl07100 

700 

CnMTlNIlE 

TR107120 

JFND  = ( JtM()*ISTA*^T)/? 

TR107110 

A(M)  = T(NiiMT*l) 

TR107130 

TRl0714n 

TO  SPACft  TO  THK  bFOTNNlNO  OF  ThF  NFXT  FUNCTION  TABDLATFO 

TH1071S0 

TRl071ft0 

linuM  r iinijK  ♦ nnphoi) 

TR107170 

1000 

CONTINUE 

TH107180 

BF  TUPn 

TH107190 

FMD 

THl072no 

SUBROUTINE  MATVeC  ( A*B.(  .N) 

TR107210 

DIMENSION  A (9) ,8 (9) ,C (9) »F (9) ,0(9) ,H(9) 

TRl 07220 

IF  (N)  10, ft. 10 

TH107230 

in 

GO  TO  (S, 6, ft, 6, ft), N 

TB107240 

ft 

DO  hi  J=1  ,9 

TR1072S0 

#<1 

F ( J) =A ( J) 

TRIO  1260 

GO  TO  70 

TR107270 

ft 

M?=l 

TR1072e0 

no  3ft  K=l,3 

TR107290 

K 1 =K,ft 

TK107300 

DO  3ft  J=K,k1 ,3 

TR 107310 

F (M2) =A ( J) 

TH1073?n 

3ft 

M?=M2,1 

TK107330 

70 

IF  (N-1)  71,71,7? 

TH107340 

71 

M4=l 

TRlU73Sn 

DO  73  J= 1 , 3 

TR1073ft0 

71 

G(J)=W( J) 

TR107370 

GO  TO  80 

TR1073H0 

7? 

M4  = 7 

TR107390 

GO  TO  (74, 74, 74, 7S, 7ft)  .N 

TR107400 

74 

DO  7ft  J=l,9 

TH107410 

7ft 

G ( J) =8 ( J) 

TR107420 

GO  TO  80 

TR107430 

7ft 

M?=l 

T81 07440 

no  6ft  K=l,3 

TR107450 

Kl=K,ft 

TR1074ftn 

DO  6ft  J=)<,Kl,.3 

TR107470 

G(M?)=R(J) 

TR1074H0 

ftft 

M?=M2, 1 

TR107490 

80 

M?=l 

THl07bOO 

no  30  M1=1.M4,3 

TR107ftl0 

no  30  K=l,3 

THlU7ft2() 

K1 =K,ft 

TR107ft30 

M3=M1 

T8l07ft40 

H (M?) =0. 

TR107ft50 

DO  20  J=K,K1,3 

THl07ftft0 

H (M2)  =H(M2)  .F  ( J)  (M3) 

TR107ft70 

?P 

M3=M3,1 

TH107S80 

30 

M?=M?, 1 

THl07ft9u 

IF  (N-1)  9ft. 9ft, 9ft 

TR1076n0 

qft 

Ml  =3 

TRl07ftl0 

GO  TO  90 

TMl07ft2(l 

Pft 

Ml  =9 

TRl 07ft30 

90 

no  91  j=i,Mi 

THl07ft40 

91 

C(J)=8(J) 

THl076ft() 

rftuwn 

TR1076h() 

FND 

TR1u767u 

SUBROUTINE  1T4B(MTAH,N,U. V) 

THl07ftrt0 

7/2/74  JOHN  F . HOLM) S 

THl07ft90 

COMMON  Y(4940) 

TP107700 

common/ TAB// (ftO) 

TR10771  (J 

COMMON/TAHI  /TAKLF  (3000) 

TR10772U 

EOUI  VAl  ENCft  ( Y (489,))  .K  r AH  ( 1 ) ) , ( Y (308  1 ) ,K  ( 1 ) ) 

•tRl  077  TO 

DIMENSION  K TAB (49) ,U ( 3) ,0 (49) 

TK10774() 

IN=KTAP  (NTAFi) 

TRl077ftO 
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KK=K ( IN) 

MFNC=1 

CAL  L FRMHAn  ( TABI.E  (KK  ) ,n.mF  nC*U.  V) 

RFTURN 

END 

SURRODTINE  FNOL  3 ( J.NN.R.L . mrr , XNF . X f Y t f' * Ot R I V . TE PM . nui RUT ) 

COMMON  YY(4940) 

DIMEN<;ION  r (3)  , Y (31  ) • YO  (31  ) ,YR  (31  ) .YC  (-<1  ) .0  (SI  ) .()M  (•’I  ,4)  t(’K  (31 .4  ) 
♦,FRPnR(31) ,YK(31 ) 

DATA  FM6»EP1 1 .M4/1 .E-h. 1 .E-1 1 *-4/ 

DATA  (C  (K)  .K  = l .■’) /?*.S.  1 ./ 

DATA  HMAX5/1.F3S/ 

NHTS=n 

FP?=n. 

N=NN 

J.I=J-? 

H = G 
HN  = H 
MK  = 1 
NRFT=M4 
JTFST=1 

IF  (JJ  .LT.  0)  jTEST=4 
IF  (XmE.EQ.  0.)  GO  TO  IS 
Fr=Y (N+3) 

FliP=ln.**(-XNF) 

EI,O=EllP«.0ni 
FM=ELn*31  .F,??776f. 
is  xr'=x 
xs=xn 
DO  20  1 = 1 
EPROR(n=0. 

20  YD ( I ) =Y ( I ) 

CALI  DEHIV (X.Y.D) 

CALI  TERM  (X*Y.n.K) 

C PPINT 

SO  CALL  OUTPUT  ( X »Y,n.FRRliH', N, I .H) 

<nnnnf**e*einnt*<nnn>*<nnnnnnnnn»**-»innt*»»»«0<nnt**tt«***«**®4nnnt************ 
* NOTE  change  made  to  ORII7TNAI  FN0L3  E OH  MOliIEI.Y 
IF (YY(?863) .FO. 1 .o)  RFTUHm 

**»'»<nnnn»******«(t*«»**o«<nnt'»*e****<nnnnj«*<nnn>**«<n»<nnnnnnn»*«********«*o 
IF  (NPFT)  FS.N0»SS 
SS  PRINT  3000,  HN 

3000  FORMATdOemEXECuriON  terminated  because  interval  of  INTEr-RATlL'N  I 
1ESS  than  l.OF  -h  TIMES  INItFPENOFNT  VAHIABLF  (X).  M =,1PF1S.7) 
STOP 

SO  RFTURM 
SS  NPR=0 

IF  (MPR  .LF.  0)  P(=Y(N*1) 

100  IF  (JTFST  .FO.  S .AND.  H .Nf.  HN ) GO  TO  4SS 
IF  (vJJ  .GE.  0)  H = HN 

IF  (Mk  .NE.  0 .OP.  JJ  .NE.  0)  GO  TO  300 

C THE  ADAMS  MOULTON  METHOD 

200  Hn24=H/24. 

JAM=0 

00  210  I=1,N 

YPI=  (SS.«Dm  ( T , I ) ♦37.»n'-i  ( I .?)  ) - (SR.*l'M(  I ,3)  ♦R.*DM  (1.4)) 

YP ( I ) =YD ( I ) ♦H0P4*YP1 
Y ( I ) =YP ( I ) 

210  CONTINUE 
X=XD*H 

CALL  OE WI V (X , Y .OM ( 1 ,4 ) ) 

DO  220  1 = 1. N 

YPI= (R.*0M ( I .4) ♦ 1Q.*Om ( I , 1 ) .dm ( T .2) ) -s."0M ( I . 3) 

YC ( I ) =YD ( I ) ♦HD24*vPI 
ERROR ( 1 ) = ( YP ( I ) -Y(  ( 1 ) ) / 1 4 . 

C THIS  ADDS  IN  A 20  COPPFCTlON 
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TR1077SS 

TR1077SO 

TP107800 

TR107810 

TR107B20 

TH107e30 

TP107840 
TR1078S0 
TH107HSO 
TH1078I0 
TR1078H0 
TR107890 
TP1U7900 
TR107910 
TRI 07920 
TRl07H3n 
TRU7940 
TRl079S() 
TH107960 
TR107970 
TH1079H0 
TR107990 
TRIOPOOO 
TRlOROlO 
TRl 0R02O 
TR10R030 
TR10H040 
THlOROSO 
TR1U80R0 
TR10H090 
TH108100 
TRIOSI 10 
TRl 0hl20 
TR10R130 
TRlUHl40 


TRIOOISO 
TRlORlftO 
TRlOPl 70 
TK10P180 
TWIOPIRO 
TR10H20O 
TR10P210 
TH10H220 
TR10P23O 
TR10P240 
TH1082S0 
TR1OP2S0 
TR10827O 
TR10n2H0 
TP108290 
TH10P300 
TR10P31fl 
TR10P32n 
TH10P330 
TR10H340 
TRl 0P3SO 
TRl OPihO 
TRl 08370 
TRl  (J83H0 
TR1(im390 
TR10H40() 
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VC(I)=YC(I)*FRROR(I) 

TR1O041O 

220 

continue 

TR108420 

IF  (XNE.NE..O)  60  TO  700 

TR108430 

GO  TO  455 

TR1O044O 

--THE  RUN6E  KlITTA  MFTHOn 

TR1O045O 

lO 

GO  TO  (301, 309, 300*309. 303) ,JTFST 

TR1O046O 

301 

DO  302  1=1, N 

TR1O047O 

YK(I)=YD(I) 

TR1O840O 

302 

CONTINUE 

THl 00490 

XOS=XO 

TH1O05OO 

GO  TO  309 

TR10H510 

303 

DO  304  1=1, N 

TR108520 

YK(I)=YC(I) 

THl 00530 

304 

continue 

TRl U054O 

xos=xn*H 

TR10R550 

30fl 

H«  = H 

TR1O056O 

H=2.*H 

TR10P570 

GO.  TO  320 

TR1O050O 

300 

X=rXO 

TH1O059O 

JAM=1 

TR1O06OO 

no  310  I=1.N 

TH10P610 

Y(i)=Yn(i) 

TH1O062O 

DK (I , n =D( T ) 

THl0rt630 

310 

CONTINUE 

TR1O064O 

IF  (JTFST  .IF.  2)  CALL  OE9 J v ( X • Y , DK ) 

THl 00650 

IF  (MK  .GT.  1 .OR.  JTEST  .OT.  1)  GO  TO  320 

TH1O066O 

DO  315  1=1. M 

TR1O067O 

DM ( I ,4) =DK ( I . 1 ) 

TR1O06HO 

3lS 

CONTINUE 

TRl 00690 

320 

DO  335  K=2.A 

TR1O07OO 

HC=H*r (K-1 ) 

TR10P710 

DO  330  1 = 1 , N 

TR10P720 

Y(I)=yD(I)  ♦ HC»OK(I«K-l) 

TR1O073O 

330 

CONTINUE 

TR10R740 

X=XO*HC 

TR1O0T5O 

call  DEKIV(X,Y,0K(1,K) ) 

TR108760 

33‘i 

continue 

TR1O077O 

Hn6=H/6. 

TR1O07BO 

DO  340  I=1.N 

TR1O079O 

YPI=DK (1,1) ♦OK ( I ,4) *2.* (UK ( I ,2) ♦OK ( I ,3) ) 

TR1O00OO 

YC  ( I ) =YD( I ) ♦HDftoYPl 

THlOHHlO 

340 

continue 

TR1O002O 

GO  TO  (360,390,370,455.370) .JTFST 

TR 100030 

300 

DO  365  I=1,M 

TR1O091O 

YP(I)=YC(I) 

TR1O092O 

36S 

continue 

TR1O093O 

JTeST=3 

TH1O094O 

GO  TO  300 

TR1O095O 

370 

DO  380  1=1. N 

TRl 00960 

Yr)(I)=YP(I) 

TR1O097O 

YP(I)=YC(I) 

TR1O890O 

300 

CONTINUE 

TR1O099O 

H=HS 

TR109000 

XOaXO+H 

TR109010 

JTFST=2 

TR109020 

IF  (NK  .EQ.  1)  GO  TO  309 

TH109030 

GO  TO  451 

TRl0904n 

390 

no  400  I=1.M 

TR109050 

ERROR  I I ) = ( YC ( I ) -YP ( I ) ) /lb. 

TR109060 

YC  ( I ) =YC ( I ) ♦ERROR ( 1 ) 

TR109070 

YP(I)=YC(I) 

TR1090H0 

400 

CONTINUE 

TR109090 

JIEST=5 

TR109100 

IF  (XNE.NE.,0)  GO  TO  700 

TR1091  10 

--ACCEPT  THF  STFP  SIZF 

TH109120 

450 

IF  (JAM  .EO.  0)  GO  TO  455 

TRl 09 130 
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IF  (MK  .EQ.  3 .AND.  JJ  .EO.  0)  GO  TO 
IF  (MK  .NE.  1)  GO  TO  303 
IF  (JTFST  .EQ.  1)  GO  in  4bb 
451  DO  452  1 = 1, N 
Y ( I ) =YD(  I ) 

45?  CONTINUE 
GO  TO  466 
45^;  DO  45<»  NQ=1  ,N 
Y0(NQ)=YC(N0) 

Y(NQ)=YD(NO) 

456  continue 

IF  (JJ  .GE.  0)  JTFST=1 

IF  (MK  .NE.  0 .OP.  JJ  .NE.  0 .OP.  NRFT  . NF . M<f)  GO  TO  466 
DO  460  l=l.N 
Dm ( I ,4)=DM (1,2) 

DM(I,?)=DM(I,3) 

DM(I,3)=DM(I,1 ) 

460  CONTINUE 
46<;  Xn  = XD*H 
466  X=XD 

IF  (MK  .EQ.  3)  MK=0 
CALL  OERIV(X,Y,n) 

DO  470  1=1, N 
DM ( I ,MK+ I ) =0 ( I ) 

470  CONTINUE 


460  FP=F 

CALL  TERM  (X,Y,0,T) 
C no  YOU  TERMINATF 


500 

IF 

(ABS(F) 

.LE.  FP6) 

GO  T 0 

800 

IF 

(FP  .EG 

. 0.)  GO  TO 

5)0 

IF 

(NRFT  . 

NF.  M4  .or. 

F«F  P 

.LT 

. FPll)  GO  TO  805 

510 

XS  = 

XD 

IF 

(MK  .NE 

. 0 .AND.  H 

.EG. 

HN) 

MK=MK+ 1 

C DO  YOU  PRINT 

600  NPR=NPR*1 

IF  (MPR  .EO.  0)  GO  TO  610 
IF  (NPR  .6F.  MPP)  60  TO  50 
60  TO  100 

610  IF  ( APIS  ( Y (N»l  ) -PC)  ,Gt.  Y(N  + 2))  GO  TO  50 
GO  TO  100 

c determining  the  step  SI/E 

700  HP  = HMAX5 

DO  760  I = 1,N 
Z=ABS(ERR0P(I) ) 

IF  (YC ( I ) .EO.O.  ) GO  TO  720 
ZZ  = YC ( I ) 

Z7=A0S (ZZ) 

IF  (EC)  720,710,705 

70S  IF  (EC  .6T.  ZZ)  ZZ=EC 
710  Z=Z/Z7 

720  IF (Z.GT.ELn.ANn.Z,LT.EliP)  GOTO  750 
HP  = 4MIN1  (HP,FM/(Z*EPn  ) ) 

GOTO760 

750  HR  = AMIN1  (HP,  1 . ) 

760  CONTINUE 

IF(HB.NE.l.)  GO  TO  765 

NHTS=n 

GO  TO  450 

765  HN=H*HP**.P 

IF  (MK  .NE.  1)  JTFST=) 

MK  = 1 

IF  (HP.LT.l.)  GOTO  770 
IF  (APS(HN)  .GT.  APS(4.*H))  = 

NHTS=0 
GOTO  450 

770  HFPS=ABS (X*FP6)  ♦ FPll 


TR109140 
TR109150 
TRl 09160 
TR1091 70 
TR109180 
TR109190 
TR109200 
TR10P210 
TR109220 
TR109230 
TR109240 
TR109250 
TR109260 
TR109270 
TR109280 
TR109290 
TR109300 
TR109310 
TR109320 
TR10R330 
TR109340 
TR109350 
TR109360 
TR109370 
TR109380 
TR109390 
TR109400 
TR109410 
TR10P420 
TR109430 
TR109440 
TR109450 
TR109460 
TR109470 
TR10O480 
TR109490 
TR1097S0 
TR109760 
TR10Q770 
TR109780 
TR109790 
TR109800 
TH109810 
TRIO  820 
TR109825 
TRl0oe30 
TR10P840 
TR10PB50 
TR10P860 
TR109870 
TRlU9b80 
TR109890 
TR10P900 

TR10Q920 
TR110130 
TRl 10140 
TRl 10150 
TR110160 
TRl 10180 
TRl 10190 
TRl 10200 
TRl 10210 
TRl 10220 
TRl 10230 
TRl 10240 
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IF  <AHS(HN)  .IT.  ABS(H/4.))  HN=H/4. 

TR110250 

IF  (ABS(HN)  .GT.  HEPS)  GO  TO  790 

TR110260 

NHTS  = NHTS  ♦ 1 

TRl 10270 

IF  (NHTS  .IE.  10)  GO  TO  78o 

TRl 10280 

NRET  = 1 

TR110290 

GO  TO  450 

TR110300 

780 

HN=SIGN(HEPS.HN) 

TRl 10310 

IF  (NHTS  .GT.  1)  GO  TO  450 

TRl 10320 

790 

IF  (NHTS  .GT.  1)  NHTS-n 

TRl 10330 

IF  (JAM  .EO.  0)  GO  TO  100 

TRl 10340 

DO  795  1 = 1. N 

TRl 10350 

YD ( I ) =YK ( I ) 

TRl 10360 

795 

CONTINUE 

TRl 10370 

XO=XDS 

TR1103H0 

JTEST=1 

TRl 10390 

GO  TO  100 

, , TR110400 

—THE  termination  LOOP 

TR110410 

800 

nret=o 

TRl 10420 

805 

IF  (NRFT  .it.  0)  GO  To  806 

TRl 10430 

H=XD-XS 

TRl 10440 

GO  TO  50 

TRl 10450 

805 

IF(F*FP.LT.0.)  GOTO  810 

TR110460 

IF  (F«FP2.L T.O. ) goto  820 

TRl 10470 

60  TO  800 

TR110480 

810 

FP2  =FP 

TRl 10490 

X 

II 

I 

TRl 10500 

GOTO  830 

TR110510 

820 

FP  =FP2 

TRl 10520 

HP  =H  + HP 

TRl 10530 

830 

nret=nret*i 

TRl 10540 

H=HP*F/ (FP-F) 

TRl 10550 

JIEST=4 

TH110560 

GOTO  300 

TR110570 

FND 

TR110580 

THIS 


TS  usFn 

limitfo 


FOP  modeling  6 PATE  GYPO  t)P  ANY  OTHEP  PROPOPTIONAl 
OUTPUT  DEVICE 


FI  = THF  INPUT 
L = The  limit  on  the  OUTPUT 
K = The  gain  of  the  device 

REAL  L*K 
KL  MT=K*EI 

IF  (AR<;(KLHT)  .GE.l  ) KLMT  = S I GN  ( L . KL  MT  ) 
END 
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DECK 

TP2 

»««»*««««« 

TR200010 

03/13/75  12. 5ft. 47  JOHN  HOLME? 

TR200012 

•«»«*««««« 

TR2000U 

OVERLAY (TRaj, 1 , 1 ) 

TR200100 

PROGRAM  SETUP 

8/2/77  JOHN  E.  HOLMFS 

TR200n0 

COMMON  Y(4940) 

TR200130 

EOUIVALENCf  (Y(2307) .NOMOO) « <Y(2312) .NMOOi 1 ) > 

TH200140 

DIMENSION  NMOrXPO) 

TR200150 

DO  1000  1=1. NOMOP 

TH200160 

L=NMOn(I) 

TR200I70 

GO  TO  (J.2.3»4»5.ft»7.8.q.l0.11.12.13.‘l4.15.1ft.l7.1B.l9.26).L 

TR200180 

1 

CALL  IMODl 

TR200190 

GO  TO  1000 

TR200200 

? 

CALL  TMOD2 

TR200210 

GO  TO  1000 

TR200220 

3 

CALL  IMOD3 

TR200230 

GO  TO  1000 

TH200240 

4 

CALL  IM0U4 

TR200250 

GO  TO  1000 

TR200260 

CALL  TMOD5 

TR200270 

GO  TO  1000 

TR200280 

ft 

CALL  IM0D6 

TR200290 

GO  TO  1000 

TM200300 

7 

CALL  rM007 

TR200310 

GO  TO  1000 

TR200320 

B 

CALL  IMOD8 

TR200330 

GO  TO  1000 

TR200340 

q 

CALI  IMOD9 

TR200350 

GO  TO  1000 

TR200360 

in 

CALL  IMODln 

TR200370 

GO  TO  1000 

TR200380 

11 

CALL  TMODll  j 

TR200390 

GO  TO  1000 

TR200400 

1? 

CALL  IMODl? 

TR200410 

60  TO  1000 

TR200420 

13 

CALL  IM0D13 

TR200430 

GO  TO  1000 

TR200440 

14 

CALL  IM0D14 

TR200450 

GO  TO  1000 

TR200460 

IB 

call  IMODIB 

TR200470 

GO  TO  1000 

TR200480 

1ft 

call  IMOOlft 

TR200490 

GO  TO  1000 

TR200500 

17 

call  IM0D17 

TR200510 

GO  TO  1000 

TR20052n 

IB 

call  TMODlB 

TR200530 

GO  TO  1000 

TR200b40 

IQ 

call  TMOOIR 

TR200550 

GO  TO  1000 

TR200560 

2n 

CA| 1 lM002n 

TR200570 

1000 

continue 

TM200580 

END 

TR200590 

TR200600 

TH200610 

SUBROUTINE  IMODl 

TB200620 

RFTURN 

TH200630 

END 

TH200b40 

TR20065n 

TR200b60 

subroutine  IMOO? 

TR200b70 
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RETURN 

TR2006eo 

END 

TR200690 

TR200700 

TR2007I0 

SUBROUTINE 

IMOD3 

TR200720 

RETURN 

TR200730 

END 

TR200740 

TR200750 

TR200760 

SUBROUTINE 

IM0D4 

TR200770 

RETURN 

TR20C780 

END 

TR200790 

TR200800 

TR200810 

SUBROUTINE 

iMons 

TR200820 

RETURN 

TR200B30 

END 

TR200840 

TR200850 

TR200880 

SUBROUTINE 

IM006 

TR200870 

RETURN 

TR200880 

End 

TR20n890 

TR200900 

TR2009IO 

SUBROUTINE 

IMOD7 

TR200920 

RETURN 

TR200930 

END 

TR200940 

TR200950 

TR200960 

SUBROUTINE 

iMone 

TR20097n 

RETURN 

TR200980 

END 

TR200990 

TR201000 

TR20I010 

SUBROUTINE 

IM0D9 

TR201020 

RE  TURN 

TR201030 

END 

TR201040 

TR201050 

TR201060 

SttBROUTINE 

IMODIO 

TR201070 

RETURN 

TR201080 

END 

TR201090 
TR201 100 
TR201110 

SUBROUTINE 

INODl 1 

TR201 120 

return 

TR201 130 

END 

TR201 140 
TH201 150 
TR201 160 

SUBROUTINE 

IMOOl? 

TR201170 

RETURN 

TR201 180 

END 

TR201 190 
TR201200 
TR201210 

SUBROUTINE 

IM0013 

TR201220 

RETURN 

TR201230 

END 

TR201240 

TR201250 

TR201260 

SUBROUTINE 

IMOO) 4 

TH201270 

return 

TR201280 

END 

TR201290 

TR201300 

TR201310 

SUBROUTINE 

IMOOIS 

TR201320 

return 

TR20I 330 
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END 

TR201340 

TR201350 

TR201360 

SUBROUTINE  IMODlft 

TR201370 

RETURN 

TR201380 

END 

TR201390 

TR201400 

TR201410 

SUBROUTINE  IMODl? 

TR201420 

RETURN 

TR201430 

END 

TR201440 

TR201450 

TR201460 

SUBROUTINE  IMODlfl 

TR201470 

RETURN 

TR201480 

END 

TR201490 

TR201500 

TR201510 

SUBROUTINE  IMODIR 

TR201520 

RETURN 

TR201530 

END 

TR201540 

TH201550 

TH201560 

SUBROUTINE  IMODPO 

TR201570 

RETURN 

TR201580 

END 

TR201S90 

TR201600 

TR201610 

OVERLAY(TRftJ,l,?) 

TR201620 

PROGRAM  INT6RT 

TR201630 

c 

7/2/74  JOHN  E.  HOLMES 

TR201640 

rOMMON  Y(4940) 

TR201650 

INTEGER  CVAR (31 ) .nVAR ( 11 ) 

TR201660 

equivalence (Y (2300) *G) , (Y(2301) .XNE) 

TR201670 

equivalence (Y (2999) .T) 

TR201680 

EQUIVALENCE (Y (2306) .ERROR) . (Y(2903) • C VAR ( 1 )).( Y ( 2934 ) ,LOCD(l) ) 

TR201690 

equivalence (Y (296S) .OVARd ) ) . (Y(2a72) .LOCC(I) ) . (Y(?P70) .NODER) 

TR201700 

EQUIVALENCE (Y (2869) .NOVAR) 

TR201710 

EQUIVALENCE (Y (2862) .TI ) , ( Y (2996) .STAGE ) 

TH201 720 

EQUIVALENCE (Y (2868) .DELT) 

TR201730 

DIMENSION  1 OCC (31 ) .LOCO (31 ) 

TB201740 

DIMENSION  r (31 ) .D (31 ) 

TR201750 

EXTERNAL  DFRIV. TERM. OUT 

TR201 760 

DO  10  1=1.31 

TR201770 

r ( I ) =0.0 

TR201 780 

10 

D(I)=0.0 

TR201790 

DO  2000  1=1 .NOVAR 

TR201800 

JJ=LOCC ( I ) 

TR201810 

JK=CVAR ( I ) 

TR201H20 

, 

C(JJ)=Y(JK) 

TR201830 

2000 

continue 

TR201840 

■ 

DO  3000  1=1. NODER 

TR20ie50 

! 

JJ=LOCD(I) 

TR201860 

JK=DVAR(I) 

TR201870 

D(JJ)=Y( JK) 

TR201880 

3000 

continue 

TR201890 

< 

IE (STAGE. GT. 0.0)  TI=T 

T = TI 

TR201900 

C(N0VAP+2)=Y(2997) 

TR201930 

C (NOVAR.3) =ERROR 

TR201940 

NN=N0VAR 

TR201950 

OELT=G 

TR20196n 

J=IEIX (Y (230?) ) 

TR201970 

L=IEIX(Y(2304) ) 

TR201980 

MPR=IEIX ( Y (2305) ) 

TR201990 

400 

CAl.L  ENOL3  (J.NN.G.  L.MpP.  XNE.  T.C.D.  DERI  V.  TERM,  (lUT) 

TH202000 

A-20 


<1- 


NSWCAVOL  TR  78^ 


END 

TH202010 

TR202020 

TR202030 

SUBROUTINE  DERIV(T»CfD) 

TR202040 

c 

7/2/74  JOHN  E.  HOI MES 

TR202050 

COMMON  Y(4940) 

TR202060 

integer  CVAR (31 ) »nVAR(31 ) 

TR202070 

EQUIVALENCE (Y(2307) .NOMOD) , (Y(2312) .NMOD(l) ) 

TR202080 

equivalence (Y (2903) ♦CVAR ( 1 )), (Y (2934) ♦LOCD( 1 ) ) 

TH202090 

EQUIVALENCE (Y (296S) ♦ DV AR ( 1 ) ) ♦ ( Y ( 2872 ) tLOCC ( 1 ) ) ♦ ( Y ( ?870 ) tNOOER) 

TR202100 

EQUIVALENCE (Y (2869) ♦ NO VAR) 

TR202110 

DIMENSION  1 OrC(31) ♦LOCO(31) 

TR202120 

DIMENSION  C (31 ) ^0(31 ) ♦NMOD(20) 

TR202130 

DO  2000  I=lfNOVAR 

TR202140 

JJ=LOCC(I) 

TR2021S0 

JK=CVAR(I) 

TR202160 

Y( JK)=C(JJ) 

TR202170 

2000 

CONTINUE 

TR202180 

DO  1000  I=ltNOMOD 

TR202190 

L=NMOn ( I ) 

TR202200 

GO  TO  (If2f3f4fs.6^7^8f9fl0flltl2fl3^l4fl5fl6fl7fl8.l9t20) fL 

TR202210 

1 

CALL  MODI 

TR202220 

60  TO  1000 

TR202230 

2 

CALL  M002 

TR202240 

GO  TO  1000 

TR202250 

3 

CALL  MOD3 

TR202260 

GO  TO  1000 

TR202270 

4 

CALL  M0D4 

TR202280 

GO  TO  1000 

TR202290 

S 

CALL  MODS 

TR202300 

GO  TO  1000 

TR202310 

h 

CALL  M0D6 

TR202320 

60  TO  1000 

TR202330 

7 

CALL  MOD7 

TR2023A0 

60  TO  1000 

TR202350 

e 

CALL  MODS 

TR202360 

GO  TO  1000 

TR202370 

9 

CALI  M0D9 

TR202380 

60  TO  1000 

TR202390 

10 

CALL  MODIO 

TH202400 

GO  TO  1000 

TR202410 

1 1 

CALL  MODll 

TR202420 

GO  TO  1000 

TR202430 

12 

CALL  MPD12 

TR202440 

GO  TO  1000 

TR202450 

13 

CALL  MOD13 

TR202460 

GO  TO  1000 

TR202470 

14 

CALL  M0014 

TR202480 

GO  TO  1000 

TR202490 

IS 

CALL  M0015 

TR202500 

GO  TO  1000 

TR202S1() 

16 

CALL  MOD16 

TR202520 

GO  TO  1000 

TR202S30 

17 

CALL  MOD17 

TR202540 

GO  TO  1000 

TH20255n 

18 

CALI  MODie 

TH202S60 

GO  TO  1000 

TH202570 

19 

CALL  MOD19 

TR202580 

GO  TO  1000 

TR202S90 

20 

CALI  MOD20 

TR202b00 

1000 

CONTINUE 

TW202610 

no  4S00  1=1 ♦NOOEP 

TR202620 

JJ=LOCO(I) 

TR202630 

JK=DVAP ( I ) 

TR202640 

D(JJ)=Y(JK) 

TR202650 

4S0n 

CONTINUE 

TH2U2660 
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RfTURN 

FNO 


SUBROUTINE  OUT ( T .C »D *EPROR » DMY » DNY X . PELT ) 

C 7/2/74  JOHN  E.  HOLMPS 

COMMON  Y(4940) 

INTEGER  CVAROl  ) ,nVAR(31)  .PLOT  (10)  .OUTNOOO) 

integer  plot 

EOUIVSLENCE ( Y (2ft«.?) .Tl ) , ( Y (2B63) .STOP) . ( Y (2864) .TS) 
E0UIVALENCF(Y(2866)  RUNN ) . (Y(2867) .PflTE) 

E(3U  I valence  (Y  (2308)  .NOOUT)  . (Y  (2332)  .NAl  ( 1 ) ) . 

♦ (Y (2392) .0UTNO( 1 ) ) 

EOU I valence (Y (231 1 ) .NOPLOT ) , ( Y (2852) .Pi OT ( 1 ) ) 

EOU I valence (Y (290  3) .CVAR ( 1 ) ) . ( Y (2934) .1 OCD ( 1 ) ) 

equivalence (Y (2965) .OVARd) ) . (Y(2872) .LOCC(l) ) . (Y{PP70) .NOPER) 
EQUIVALENCE  (Y  (3044)  .TITLd)  ) 

E0UIV4LENCE(Y(2869) .NOVAR) 

EQUIVALENCE  (Y  (3043)  .PLCT)  . (Y(2362)  .HOl  d ) ) . ( Y ( 2377  ) .H02d)  ) 
DIMENSION  I OCC (31 ) .LOCn (31 ) 

DIMENSION  TITL  (7)  .HDl  (15)  .H02d5)  .EORl  (31 ) .E0p2(31  ) 

DIMENSION  NAl (30) .C(3l) .D(31 ) 

REAL  PLdO) 

1 EPPMAT (3X.4HRUN  . F 12.2 . 1 OX . 7 A 1 0 .6X . 7HE ORMAT  .11) 

2 E0RMAT(3X.F12.2.9?X.5hdaGE  ,13) 

4 FPRMATdHl) 

5 FORMAT dH0,2X.5HTIME  ,15A8/) 

**  K?  = NO.  LINES  STORED  FOR  CURRENT  PAGE 
**  NUM  = NO.  formats 
♦*  NP  = NO.  PAGES  OF  EACH  FORMAT 
**  J.J2  = SUBSCRIPTS  FOR  OUTPUT  ARRAYS 
DO  3000  I=l,NODER 
JJ=LOCD(I) 

JK=DVAR ( I ) 

Y ( JK) =D ( JJ) 

3000  CONTINUE 

DO  2000  1=1. NOVAR 
JJ=LOCC(I) 

JK=CVAP ( I ) 

Y(JK)=C(JJ) 

2000  CONTINUE 

Y(2868)=DEl  T 
Y (2999) =T 

IF(T-TI)  105,100.105 
100  CONTINUE 
NP=0 
K7  = 0 
J2=0 
J=0 

105  CALL  PROCESS 

»»»  «;tOPE  output  IN  Y(3100+)  ARRAY  UNTIL  FULL  PAGE  IS  ACCUMULATEP 
IIS  K7=KZ+1 
200  J=J»1 
J2=J2.1 
Y(J.3100)=T 
Y(J2*3980)=T 
DO  330  1=1. NOOUT 
IN=OUTNO(I) 

IF(I.IF.15)  J=J.l 
IF(I.LE.15)  Y(J*3100)=Y(IN) 

IF(I.GT.15)  J2=J2*1 
IF(I.GT.15)  Y(J2*I980)=Y(1N) 

330  CONTINUE 

599  IF  (NOPLOT. FQ.O)  GO  TO  ■<31 
DO  600  1=1. NOPLOT 
1I=PL0T(I) 


TR202670 

TR202680 

TR202690 

TR202700 

TR202710 

TH202720 

TR202730 

TR202740 

TH202745 

TR202750 

TR202770 

TR202780 

TR202790 

TR202800 

TR202810 

TR202820 

TR202830 

TR202835 

TR202840 

TR202850 

TR202860 

TR202870 

TR202880 

TR202890 

TB202910 

TR202920 


TR202940 

TR202950 

TR202960 

TR202970 

TR202980 

TR202990 

TR203000 

TR203010 

TR203020 

TR203030 

TR203040 

TR203060 

TR203080 

TR203100 

TR203110 

TH203120 
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niMFNSIONFD  FOR  5i>  LINFS 


GO  TO  IR 


600  Pi  <I)=r<  II) 

PLCT=PLCT*1 

WRITE  ( 19)  RUIJN.T.  (PL  ( I ) 1 1 = 1 .NOPLOT  ) 

331  IF(STOP)  33?, 33?. 15 
33?  IF(KZ.LT.5n  GO  TO  333 
***  WRITE  PAGE  OF  OUTPUT  IN  EACH  FORMAT . .ARRAYS 
15  NP=NP.l 

IF (NP.GT. 1 ) GO  TO  ?? 

NUM=(lA.N0OUT)/15 
FORI (1 )=9h(lX,F7.?, 

FOR?(l )=9h(lX,F7.?, 

DO  19  N=l,15 
I = ?*N 

FORI (I)=HDl (N) 

FORI (I*1)=1H, 

IF (NOOUT.GF.N) 

FORI (I)=10h 
FORI (I.1)=1H 
IR  FOR?(I)=HD?(N) 

F0R?(I*1)=1H, 

IF ( (NOOUT-15) .GE.N) 

FOR?  1 1 ) =10h 
F0R2(1*1)=1H 
19  CONTINUE 

IF(NOOUT.GF.IS) 

IF (N00UT.lt. 15) 

IF (NOOUT.LF. 15) 

IF  (NOOUT.GT. 15) 

??  DO  30  L=1.MUM 
WPITE(6.1) 

WRITE (6,?) 

IF  (L.GT.l) 

?0  WPITE(6,5) 

M=J.3100 
WRITE (6.F0P1 ) 

GO  TO  ?6 

?1  WRITE(6,5)  (NAl(l) 

M=J?*3980 

WRITE (6, FOP?)  ( Y ( I ) » I=T981 .M) 

?6  WPITE(6,4) 

30  CONTINUE 
J = 0 
J?  = 0 
K7  = 0 

GO  TO  333 
333  RETURN 
END 


SUBROUTINE  TERM(T,C.n,F) 

7/?/74  JOHN  F.  HOLMES 
COMMON  Y(4940) 

integer  CVAP(31) ,0VAR(31) .STPNO(IO) 

EOUIVALENCF  (Y(?Rh?).TT) 

F(TUIVALENCF  (Y  (?853)  .STOP) 

EOUI VAl ENCF ( Y (?31 0) .NOSTOP) , (Y (?8??) .STPNO ( 1 ) ) . ( Y ( ?P3? ) .SUP ( 1 ) ) . 
♦ (Y (?84?) ,Sl  0( 1)  ) 

EIjUIVALENCF  ( Y (?9  0 3)  .CVARd  ) ) , (Y(?934)  .(  OCU  ( 1 ) ) . (Y(P»t>9l  .NO  VAR) 
FOUI VALENCF (Y (?965) . DV AR ( 1 ) ) . ( V ( ?e 7? ) . I OCC ( 1 ) ) , ( Y ( ?P70 ) .NOOER) 
EOUI valence (Y (?009) ,LRP ( 1 ) ) 

DIMENSION  I OCr ( 31 ) .LOCO (31 ) 

DIMENSION  SUP ( 10) ,SLO( 10) .r (31 ) ,0 (3) ) ,M (9) .RPT (9) ,RPl (9) 

PFAL  LRP(9) 

Y (?999) =T 
F=1.0 

DFLA=LRP  (1 ) *1  MP  (5)  *LPP  (9)  ♦(  RP  (?)  *LRP  (6)  *LPP  ( 7)  ♦LRP  ( ■<)  * 


GO  TO  19 


FORI (31 ) =1H) 

FORI (NOOUT*?*! ) =1H) 
FOP?(l )=9H(1X) 

F0R2( (NOOUT-15) *?*1)=1H) 


RUNN. (TITL (I) ,1  = 1 .7)  ,L 

DATE.NP 

GO  TO  ?) 

(MAI ( I ) , 1=1 . 15) 


(Y ( I ). 1=3101 .M) 


.1=14.30) 


TR?03230 

TH?03?40 

TR20324? 

TR203243 


TR203260 


TH203290 

TR203320 

TR203330 

TR203340 

TR203370 

TR?03530 

TR203540 

TR203550 

TR203580 

TR203590 

TR204100 
TH2041 10 
TH?04l?0 
TR?04130 
TR?04140 
TR204150 
TR?04lbn 
TH20A170 
TR?04l«)) 
TR?04190 
TH?04?00 
TW?U4?10 
TH204??(I 
TR?0423n 
TR?04?40 
TR204250 
TR?04?b0 
TP204?7n 
TH?04?8H 
TK204?90 
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♦ LRP(B)  *LHP  (♦)  -LRP(7)*LBP  (b)  *LRP(3)-LP^^  (»)  *LPP(f>)  *L‘»P<  1 ) 
♦LPP(9)*LRP(?)*LRP(4) 

IF(DElA)  50«60.SO 
50  CONTINUE 

C***»TO  INSURE  orthogonality  OF  LRP  MATRIX  ***** 

DO  9 I=lt3 

CALL  MATINV(LRP»RPTtB) 

CALL  MATINV<RPT.RtRPI) 

DO  8 J=l*9 

LRP( J) = (LRP( J) ♦RPT ( J) ) /?. 

CONTINUE 
60  CONTINUE 

DO  15  I=1.N0DER 
JJ=LOCO(I) 

JK=DVaR ( I ) 

Y(JK)=D(JJ) 

15  CONTINUE 

DO  10  I=ltNOVAR 
JJ=LOCC<I) 

JK=CVAR(I) 

Y(JK)=C(JJ) 

10  CONTINUE 

IPRVA=IFIX (Y (?99B) ) 

C(NOVAR*l)=Y(IPRVA) 

SGN=1 . 

DO  20  1=1. NOSTOP 

IN=STPN0 ( I ) 

HALF  = ABS(SiiP(I)-SLO<I)  )/2. 

SM1D= (SUP ( I ) ♦SLO ( I ) )/2. 

FF=HAIF-ABS ( Y ( IN) -SMI  0) 

IF(FF.LT.O.)  SGN=-1. 

IF(FF.LE.O.)  STOP=l. 

F=F*ABS(FF) 

20  CONTINUE 
F=F*SGN 
3?  RETURN 
END 


SUBROUTINE  MODI 

RETURN 

END 


SUBROUTINE  MOn2 

RFTURN 

END 


SUBROUTINE  M003 

RFTURN 

END 


SUBROUTINE  M004 

RFTURN 

END 


SUBROUTINE  MODS 

RFTURN 

END 


SUBROUTINE  MOD6 
RFTURN 

A-24 


TR20A300 

TR204310 

TR204320 

TR204330 

TR204340 

TR2043b0 

TR204360 

TR204370 

TR204380 

TH204390 

TR204400 

TR204410 

TR204420 

TH204430 

TR204440 

TR2U4450 

TR204460 

TR204470 

TR204480 

TR204A90 

TR204500 

TB204510 


TH204520 

TR204530 


TR204590 

TR204600 

TR204610 

TR204620 

TR204630 

TR204640 

TR204650 

TR204660 

TR204670 

TR2046flO 

TR204690 

TR204700 

TH204710 

TR204720 

TH204730 

TH2047*0 

TR204750 

TR204760 

TR204770 

TR204780 

TH204790 

TR204B00 

TR204bl0 

TH204B20 

TR204830 

TR204b40 

TR204850 

TR2048N0 

TH204870 

TR204880 

TR?0<.890 
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EMO 


SUBROUTINE  MOOT 

RETURN 

END 


SUBROUTINE  MODB 

RETURN 

END 


SUBROUTINE  MOOR 

RETURN 

END 


SUBROUTINE  MOniO 

RETURN 

END 


SUBROUTINE  MODI  I 

RETURN 

END 


SUBROUTINE  M0012 

RETURN 

END 


i 


SUBROUTINE  MODI  3 

RE  TURN 

END 


SUBROUTINE  M0DI4 

RETURN 

END 


SUBROUTINE  MOniS 

RETURN 

END 


SUBROUTINE  M0016 

RETURN 

ENO 


SUBROUTINE  MOnl7 

RETURN 

END 


SUBROUTINE  MODlfl 

RETURN 

END 


SUBROUTINE  MODIR 

return 

END  A-26 


TR204900 

TR204910 

TR204920 

TR204930 

TR204940 

TW204950 

TR204960 

TR204970 

TR2049tiO 

TH204990 

TR20N0O0 

TR20S0I0 

TR20S020 

TR20S030 

TR20S040 

TR20S050 

TR20S060 

TR20S070 

TH20S080 

TR20S090 

TR20B100 

TR20S110 

TR20S120 

TR20S130 

TR20B140 

TR20B1S0 

TR20B160 

TH20B170 

TR20^I80 

TP20SI90 

TH20S200 

TR20B210 

TR20S220 

TR20‘^230 

TR20S240 

TR205260 

Tft20S260 

TR20B270 

TN20B280 

TH20S290 

TH20S300 

TR20B3I0 

TR20S320 

TR20S330 

TH20S340 

TR20S350 

TR20S3B0 

TR20E370 

TR2OS3H0 

TR20S390 

TH205400 

TH20S4I0 

TR205420 

TR20S430 

TR20B440 

TR20S450 

TH20E4f.0 

TR20B470 

TR2054H0 

TH2UB490 

TH20Bbno 

TREUbblO 

TR20SS20 

TR20''SJO 

TH20bS4O 

TK2(IB550 
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TH20SS60 

TR205S70 

SUBROUTINE  M0020 

TR20SSH0 

RETURN 

TR20t590 

END 

TR20S600 

TR20S610 

TR20SG20 

SUBROUTINE  PROCESS 

TR20S630 

RETURN  ' 

TR20SG40 

END 

TR20S6S0 

TR20S660 

OVERLaY(TRaj,?.0) 

TR20S670 

PROGRAM  TRJPITS 

TR20S680 

c 

7/2/74  JOHN  F.  HOLMFS 

TR20M690 

COMMON  Y(494n) 

TR20S700 

CALL  MYPLOT 

TR20K710 

FND 

TR20S720 

SUBROUTINE  MYPLOT 

TR20S730 

c 

n/13/74  JOHN  f.  HOI  mFS 

TR20S740 

COMMON  Y(4940) 

TR20S76() 

INTEGER  PLOT(lO) .PLOT 

TR20S760 

EOUI VALENCF ( Y (2311 ) t NOPLOT) , (Y (2B52> .Pi  OT ( I ) ) 

TR20'.770 

EOUIVALENCF (Y (3043) .PLCT) 

TR20S77S 

RFAL  PL ( 10) 

TH20S780 

REWIND  19 

TR20S78S 

c 

TR20S790 

c 

THE  INFORMATION  TO  BE  PLOTTED  SHOULD  HE  READ  AS  FO|  1 OWS. 

TR20S800 

c 

TR20S810 

c 

DO  20  J=1.PLCT 

TP20S81S 

c 

RF AD ( 19)  RUNN.T . (PL ( I ) • 1 = 1 .NOPLOT ) 

TR20SH20 

c 

2n 

CONTINUE 

TR20S82'j 

c 

2n 

CONTINUE 

TR20S830 

c 

PI CT  ARE  The  number  OF  DATA  POINTS  TO  ME  PLOTTED 

TB20S83S 

c 

SUPPLY  YOUP  OWN  GOULD  CALLS 

TR20SH40 

c 

TR20S8S0 

RPTURN 

TR20MeG0 

FND 

TH20S87fl 
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■> 

APPENDIX  B 

u 

FIXED  STORAGE  ASSIGNMENTS 
Y ARRAY  LOCATIONS  2300  - 4940 


NOTE:  I 

I 

An  asterisk  after  the  Y array  location  of  the  parameter  Indicates  I 

that  the  parameter  Is  to  be  placed  In  the  data  as  an  Input  parameter.  i 


B-1 


A ' 
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PARAMETER 

IDENTIFICATION  LOCATION 

IN  Y ARRAY 

G 

2300 

* 

XNE 

1 INTEGRATION 

2301 

* 

J 

[ CONTROLS 

2302 

* 

NN 

) SEE  PAGE 

2303 

* 

L 

2304 

MPR 

1 

2305 

* 

ERROR 

2306 

* 

NOMOD 

Number  of  modules 

2307 

NOOUT 

Number  of  output  parameters 

2308 

NOIN 

Number  of  code  3 parameters 

2309 

NOSTOP 

Number  of  STOP  conditions 

2310 

NOPLOT 

Number  of  plot  variables 

2311 

NMOD(l) 

2312 

+ 

MOD  controls,  see  subroutine 

4 

INPUT 

NMOD(20) 

2331 

RNMEl(l) 

1 

2332 

Header  controls,  see  subroutine  4 

, INPUT 

RNME1(30) 

1 

2361 

HDl(l) 

1 

2362 

+ 

1 

4 

HD1(15) 

Header  controls,  see 

2376 

subroutine  INPUT 

HD2(1) 

2377 

+ 

4 

HD2(15) 

2391 

OUTND(l) 

2392 

Output  controls,  see  sub- 

4 

routine  INPUT 

OUTNO(30) 

2421 

INNO(l) 

1 

2422 

+ 

, Input  controls,  see  sub- 

4 

1 routine  INPUT 

INNO(200) 

1 

2621 

VALVE(l) 

1 

2622 

4- 

Input  controls,  see  sub- 

4 

1 routine  INPUT 

VALVE(200) 

1 

2821 

STPNO(l)  J 

1 

2822 

1 

) Stop  (termination)  controls. 

4 

STPNO(IO)  j 

1 see  subroutine  INPUT 

2831 

SUP(l) 

1 

2832 

4 

[ Upper  bound  termination 

4 

1 locations 

SUP (10) 

2841 

SLO(l) 

1 Lower  bound  termination 

2842 

( locations 

4 

SLO(IO) 

2851 

B-2 
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PARAMETER 

IDENTIFICATION 

LOCATION  IN  Y ARRAY 

PLOT(l)  } 

2852 

+ t 

Plot  variable  locations 

PLOT (10)  ; 

2861 

TI 

Initial  time 

2862 

* 

STOP 

Termination  parameters 

2863 

TS 

286A 

PRFR 

Print  control,  see 

2865 

* 

RUNN 

RUN  Number  ) Printed  on  each 

2866 

* 

DATE 

Date  f page  of  printout 

2867 

* 

DELT 

Integration  time  step  (see) 

2868 

NOVAR 

Number  of  dependent  variables 

2869 

NODER 

Number  of  derivatives 

2870 

NOTAB 

Number  of  tables 

2871 

LOCC(l) 

2872 

Locations  of  dependent  variables 

+ 

LOGO (31) 

In  C array 

2902 

CVAR(l)  ) 

1 Y array  locations  of 

2903 

^ dependent  variables 

+ 

CVAR(31)  j 

2933 

LOCD(l) 

1 Locations  of  derivatives 

2934 

I 

> in  D array 

+ 

LOCD(31)  ' 

\ 

2964 

DVAR(l)  ] 

1 Y array  locations  of 

2965 

+ ' 

^ derivatives 

4- 

DVAR(31)  ) 

1 

2995 

KEND 

Size  of  table  array 

2996 

— 

Print  increment  when  MPR=0 

2997 

* 

PRVA 

Print  variable  when  MPR*0 

2998 

* 

T 

Time  (sec) 

2999 

* 

RE 

Earth's  radius,  default  value  Is 

3000 

* 

20925631  ft 

WE 

Earth's  spin  rate,  default  value 

3001 

* 

is  7.29211508  x 10“^  rad/sec 

3002 

3003 

D 

Aerodynamic  reference  length  (ft) 

3004 

* 

A 

^ Aerodynamic  reference  area  (ft^) 

3005 

* 

DXG 

t 

3006 

* 

DYG 

j Center  of  gravity  location  w/r 

3007 

* 

DZG 

to  geometric  axes,  see  Figure  7 

3008 

* 

IXX 

3009 

* 

lYY 

Principal  moments  of  inertia 

3010 

* 

(slug-ft2)  ^ 

IZZ 

3011 

* 

MS 

Vehicle  mass  (slug) 

3012 

* 

H 

Vehicle  altitude  (ft) 

3013 

* 

TAUR 

Longitude  (deg) 

3014 

* 

PSIR 

. Latitude  (deg) 

3015 

* 

^XY 

^XZ 

3016 

* 

( Cross  products  of  Inertia 

3017 

* 

IyZ 

B-3 

3018 

* 
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PARAMETER 


PLCT 

TITL(l) 

+ 

TITL(7) 

K(l) 

+ 

K(49) 


KTAB(l) 

+ 

KTAB(49) 


IDENTIFICATION  LOCATION  IN  Y ARRAY 

3019 

3020 

3021 

3022 

3023 

3024 

3025 

3026 

3027 

3028 

3029 

3030 

3031 

3032 

3033 

3034 

3035 

3036 

3037 

3038 

3039 

3040 

3041 

3042 


Number  of  data  points  being  plotted  3043 

Storage  location  for  3044 

title  + 

3050 

Intermediate  table  3051 

Identification;  see  subroutine  ITAB  + 

3099 

3100 

Output  array  storage  + 

4889 

Assigned  array  locations  4890 

for  tabulated  functions  + 

4938 

4939 

4940 
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APPENDIX  C 

IMODl;  INITIAL  DIRECTION  COSINE  MATRIX  FOR 
3DOF  OVER  A ROTATING  SPHERICAL  EARTH 


As  mentioned  previously,  MODI  contains  the  direction  cosine  matrix. 
Since  it  would  be  unwieldly  to  read  in  the  initial  values  of  the  direction 
cosine  matrix,  IMODl  is  used  to  calculate  these  values  from  the  initial 
latitude  and  longitude  which  must  be  read  from  the  input  cards.  This 
transform  matrix  for  transforming  vectors  from  the  inertial  axes  to  the 
local  axes  is  presented  as: 


COSl^J^ 

0 

sinij/j^ 


siniiiRsinXj^ 

COSTp 

-cosii'j^sini 


-sinii'^cosi 

sintj^ 

COSlJ/pCOSTp 


Input 


Output 


IMODl  Parameters 


Parameter 

Units 

Location  in  Y Array 

^’r 

Degree 

3015 

^R 

Degree 

3014 

r- 

2000-2008 

C-1 
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c 

c 

SUBROUTINE  IMOOl 

TR300100 

3/10/78 

c 

TB300110 

c 

This  routine  calculates  the  initial  direction  cosine 

MATRICE  LRI 

TR300150 

c 

FOR  3nOF  SIMULATIONS 

TR3U0160 

c 

TH300170 

c 

IT  REQUIRES  THE  FOLLOWING  INITIAL 

CONDITIONS  ON  CODE 

3 CONTROL 

TR3001fl0 

c 

CARDS 

TR300190 

c 

TR300200 

c 

TAUR  = Y(30l4)  , LONGimOE  (DEG) 

c 

PSIR  = Y (301':)  . LATITUDE  (0E6) 

c 

COMMON  Y(4Q40) 

TR3U0230 

equivalence (Y(301S)»PSTR)»(Y(30I4) 
EOUIVAl.ENCF  (Y  (200(1)  «LRI  (D) 

,TAUR) 

TR3002ftO 

REAL  LRL(9i 

TH300270 

CALL  SENC0S(PSIR.SP*CP,0) 

TH300280 

CALL  SENCOS (TAUR.ST*CT.O) 

TR300?90 

LRL ( 1 ) =CP 

TR300300 

LHL (?) =0. 

TR300310 

LPL (3) =SP 

TB300320 

LPL (4) =SP*ST 

TR300330 

LRL  (S)=CT 

TR300340 

LRL (ft) =-CP*ST 

TR300350 

LRL (7) =-SP*CT 

TH300360 

LRL (fl)=ST 

TR300370 

LRL (9) =CP«rT 

TR3OO30O 

RETURN 

TB300390 

END 

TR300400 

TR3004IO 

TR300420 

TR300420 
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APPENDIX  D 

MODI;  DIRECTION  COSINE  MATRIX  FOR  3DOF  OVER  A 
ROTATING  SPHERICAL  EARTH 


The  purpose  of  this  module  is  to  calculate  the  direction  cosine  matrix  for 
a vehicle  flying  a 3D0F  particle  trajectory  over  a rotating  spherical  earth. 

The  longitude  and  latitude  of  the  vehicle  are  calculated  from  the  inertial 
coordinates  as. 


The  direction  cosine  matrix  which  allows  for  transforming  vectors  from  the 
Inertial  axes  to  the  local  axes  is  calculated  as. 


t^RLl  = 


COSlJij, 

0 

siniiir 


sinil/^sinij^ 


COST, 


-cosij^j^sinTj^ 


-sinij^j^cosTj^ 


sint. 


COS\|ipCOSTj^ 


The  transfer  matrix  [S-^rI  expressed  as 


LRJ 


‘'RL' 


Input 


Output 


Parameters 


X 

Y 

Z 


R 

R 

R 


R 

^'r 

t*^RLl 

t^LRl 


MODI  Parameters 


Units 

Location  in  Y Array 

feet 

803 

feet 

804 

feet 

805 

radian 

3014 

radian 

3015 

2000-2008 

2039-2047 

D-1 


non  N-  n o o 
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C 

c 

c 

c 

c 

c 

c 


0 


SIIPROIITINE  MODI 


TR300930 

TR300940 


THIS  ROUTiMt  CALCULATES  THF.  DIRECTION  COSINF  mATRICF  LRL  FOR  A TR3U09AO 
3DOF  SIMULATION  OVFR  A ROTATINR  SRHFPIfAL  FARTH  TR300y90 

TR301000 

3/10/78 


COMMON  Y(<»940) 

FOUIVALENCF ( Y (803) .XR) , ( Y (804) « YP ) . ( Y ( F 05 ) t ZR ) TR301020 

FOUI VALENCF (Y (3015) »HS I R ) » ( Y ( 30 1 4 ) .TAUR) 

FOUIVALENCF ( Y(2000) fLPI  (1))  TR301040 

FOU I valence (Y (?0 39) ♦LLR(l) ) TR301050 

RFAL  LRL(9),lLR(9)  TR301060 

DIMENSION  R(9)  TR30I070 

LATITUDE  And  lonoitude  calculated 

CALL  ARKTAN(-YR.ZHtTAUR.l ) TR301080 

CALL  ARKTAn(XR.  (SrjRT(ZW**2  + YR**2)  ) .PSlRtl)  TR301090 


INERTIAL  TO  LOCAL  AXES  TRANSFER  MATRIX 


CALL  SFNC0S(TAIJR.STAR.CTAR,1) 

TR30I 100 

CALL  SFNCOS (PSIR.SPSR.CPSR, 1 ) 

TR301110 

LPL ( 1 ) =CPSp 

TM301 120 

LRL(2)=0. 

TR301130 

LRL (3) =SPSR 

TR3U1140 

LPL (4) =SPSP«SrAR 

TR301150 

LPL  (5) =CTAP 

TR301160 

LPl (6) =-CPSP*STAR 

TR301 170 

LPL (7) =-SPSR*CT4R 

TR301 180 

LPL (8) =STAP 

TR301 190 

LPL (9) =CPSR*CTAR 

TR301200 

CALL  mATINV (LRL.H.LLR) 

TR301210 

RF  TURN 

TR3U1220 

FnD 

TR301230 

TR301240 

TR301250 

TR301 260 

D-2 
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APPENDIX  E 

IMOD2;  INITIAL  CONDITIONS  FOR  A 3DOF  OR  6DOF  TRAJECTORY 


The  purpose  of  this  module  is  to  calculate  the  initial  values  for  the 
inertial  coordinates  and  velocities  of  the  vehicle  flying  over  a spherical 
earth. 

Since  the  initial  altitude  latitude,  and  longitude  of  the  vehicle  are 
specified  in  the  input  conditions,  the  inertial  coordinates  can  be  expressed  as, 

= sinit/j^(h  + Rg) 

= -cosii<j^(h  + Rj,)sinTj^ 

Zg  = cosi()R(h  + Rg)cosTj^ 

Likewise,  the  velocity  components  with  respect  to  the  inertial  axes  can  be 
expressed  as  functions  of  the  initial  inputed  flight  path  angles  and  resultant 
velocity  as. 


VgCOSYgCOSYj^ 

^R 

= 

-VgCosYgSinYj^  - Wg(Rj,  + h)cosi|jj^ 

. zr. 

The  initial  starting  point  of  the  vehicle's  flight  path  as  projected  onto  the 
earth's  surface  is  defined  as. 


This  IMOD  serves  for  M0D3  and  M0D9  also. 


E-1 
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IMOD2  Parameters 


Parameter 


Units 


Location  in  Y Array 


h 

feet 

3013 

Re 

feet 

3000 

Ve 

(Rg  default  value  = 20925631  feet) 
feet/second 

610 

degree 

2208 

"^E 

degree 

2209 

radians/sec 

3001 

(tOg  default  value  = 0.0000729211508 

rad /sec) 

"r 

degree 

3014 

^R 

degree 

3015 

""r 

feet 

803 

Yr 

feet 

804 

Zr 

feet 

805 

^R 

feet/second 

800 

Yr 

feet /second 

801 

Zr 

feet/second 

802 

R 

feet 

818 

^^i 

R,„ 

feet 

819 

o r>  .-)  V r>  n 
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r 

c 

r 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


r 

r 

r 


SliBrtOllTlNE  lMni1? 

3/10/7H 

This  boutimf  raimiATFs  ihf  imtial  vai  uts  of  xBtVP.ZB 

AND  XPD.  YPD.  and  7HI).  AND  THt  INITIAL  EAhTH  FOSIFI'^NS 
RTF  I.  and  F>PFI 


ThF  FDlLOxlNr-  INIIIAl  rONl'ITIONS  AWT  WFUUIRFD  ON  f one  i 
CONTROL  CACPS 


PSIR  = Y(301B)  , I ATITliOF  (utO) 

TA()R  = Y (301<. ) . I ONOITUDF  (DFI’F 

H=Y(3nl3)  , ALTTTtJDtFFr) 

RF=Y(3000)  . RAOIUS  '1F  FARTHFFT),  DtFAUlT  VAUJE 

WF=Y(300I)  . EAkTHS  spin  kA TF ( PAo/ Sf C ) . DFFADLT 

0.0n0072a?l ISOrt 

VF=Y(M0)  . TNTTFAL  YFLOCITY  W/R  to  FARTh  (fPS) 
GAMAE  = Y (2204)  . velocity  FIEVATION  ANGLE  (OF.O) 

GAMAH  = Y (2208)  . VELOCITY  HEAlTING  ANOl  E (DEG) 


= 20425F3I. 
WAIUE  = 


TR3(<0430 

TR3UU44t) 

TR300*Ei0 

TR3004RO 

THaOObOO 

TR3O0bI0 


COMMON  Y(4040) 

EOl'I  VALENCF  (Y(Rm).XR).(Y(8l)4).YB),(Y(80i>).7R) 
Fcnl  VALENCf  ( V (?n  TO)  .(  L R ( I > ) • ( Y ( 20('0  ) •(  p(  ( I ) ) 
EDO  I VALENCF  ( Y ( TO  IS)  •PS’R  ) • ( Y ( ill)  4 ) . TfldP  ) 

F 00 1 VALENCF  ( Y ( Tnoi  ) *K'E  ) 

pool VALENCF  ( Y (Ml  8)  if  1 ) . ( Y (E 10)  .RPF 1 ) 

FOOT VALENCF (Y(3nnO),RF),{Y(30l3).8) 

Eooi  Vai.EnCf  (Y  (221  ’)  »0T  ) , (Y  (2214)  ,0P) 

FOOl VAI  ENCF ( Y (6Sh)  .Ml ) , ( Y ( TOl 2)  .MS) 

EFJO  I VALENCF  (Y  (??0H)  .GA'lAH)  . (Y  (220S)  .(-.AmAE  ) 

PFAL  LLH  (9>  .1  PI  (R)  .LR 
DIMENSION  P(Q) 

ENTRY  IMOIV. 

FF'TPY  JMOOo 

IMTlAt  INFPTIAL  POSIITON 

call  SFNCOS (TAOW.STAR.cTAR.U) 

CALL  SFNCOS (PSTW.SPSR.rPSR.d) 

XP  = SPSP«  (Fi.PF  ) 

LS=CPSR* (H.RF I 
YP=-( poSTAP 
ZP=LR*CT AR 

INITlAl  INFWTTAI  VELOl.TTIEs 

call  SFNCOS  (f-AMAF.SE  .('f, (I) 

CA[  L SFNCOS  (OAMAH.SN.C  s.  (j ) 

Y (512) =Y (61 0) *CF®CM 

Y (SI  3)  =-Y  ( A1  0 ) *CF  <>SH  « 

Y (514) =Y (hlOl «SF 

Y (SI  3)  = Y (51  3 1 - WF«  (PE  .F- ) »(>SP 
CALL  MAT  IN'/  ((  P|,  ( 1 ) .M.l,  I P ) 

CALL  mat  VEc  (I  I P.Y  (S12)  .Y  (MiXil  .8) 

Do  20  T=1.T 

Y ( T *51  'E)  =Y  ( 1 fSI  n 

INITIAI  RAfK.FS 

PAD  = T.  1415‘)?ASTS8g/l8(i. 


TRSUOSSO 
TR30nS  h) 

TR3(/f)5<>i> 

TR30n600 

TR30061O 

TH3(/062(y 

TR300630 

TP30064(t 

TR300650 

TF(3u066(l 

TR3on66S 


TR300670 

TR3un68n 

TR30069n 

TP300700 

TW3007in 

TH30072() 


TW3on73() 

TH30n74t1 

TR3on7SO 

TR30076il 

Th30P77() 

Tp.TonfHil 

Ik3007'Ti) 

TRTuoMno 

TR30ft8h() 

Th3o087i( 


TRTOPlFlO 


E-3 
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APPENDIX  F 

MOD2;  3D0F  PARTICLE  TRAJECTORY  ALONG  A 
PROGRAMMED  FLIGHT  PATH 


The  purpose  of  this  module  is  to  calculate  the  transformation  matrix  for 
a vehicle  flying  along  a 3D0F  particle  trajectory  over  a spherical  earth  where 
the  elevation  angle  of  the  velocity  vector  is  preprogrammed.  This  elevation 
angle,  is  taken  from  a table  of  Yj,  as  a function  of  time  which  is  to  be 

supplied  by  the  user  as  Table  Array  Number  10.  The  heading  angle,  Yjj»  is 
maintained  constant.  The  altitude  is  calculated  as, 

/ 2 2 2 
h = +Y^  -R, 

The  local  velocities  with  respect  to  the  earth  are. calculated  as. 


LE 

L “ 


The  resultant  velocity  is  then  calculated  as. 


2 2 2 
/ +V  +V 

\e  \e  ^LE 


The  matrix,  necessary  for  transforming  vectors  from  thfe  principal  axes 

to  the  local  axes  is  written  in  terms  of  and  y_  as. 


COSYj^COSYg 


-sinYj^cosYg 


-cosYjjSinYj. 


sinYj^sinYg 
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Input 


Output 


MOD2  Parameters 


Parameter 

Units 

Location  in  Y Array 

degree 

2209 

degree 

2208 

feet 

803 

feet 

804 

2r 

feet 

805 

Xr 

feet/second 

800 

Yr 

feet /second 

801 

"r 

feet /second 

802 

feet 

3000 

2057-2065 

h 

feet 

3013 

le 

feet/second 

512 

feet /second 

513 

^LE 

v„ 

feet/second 

514 

^LE 

V. 

feet/second 

610 

I 


c. 


! 

L 


t 


F-2 


r 


f. 


i 


r>ODnr>o  non  non 
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SUBHOIJTINE  MOD? 

C 

C pOSITTOM  MOntiLt  EOR  a 3D0F  PARTICIE  TRAJ.  WITw  PROCtRAMMED 

C PITCH  ANGLE 

C 

C 3/10/78 

C 

C GAMAF  (DEG)  IS  TAHULATFO  AS  A FUNCTION  OF  T I ML  (SEC)  IN 

C TARLE  AorAY  no.  10 

C GAMAH  (DEG)  IS  A CONSTANT  IN  LOCATION  Y(2^0fl) 

C 

COMMON  Y(404n) 

FOUIVALENCF  (Y(3ni3).H) 

Einil valence (Y (803) .XR) , ( Y (804) .YR) . ( Y (805) .ZR) 

EOUIVALENCF ( Y ( 301R) *PSrR) » ( Y(3014) .TAUP) 

EO))I  VALENCE  (Y(3n01)«V)t)«iY(3000)*RE) 

EOUIVALENCF (V(S12),VXLE)«(Y(513)tVYlF).(Y(514),VZLF) 
EOUIVALENCF  ( Y (POPS))  »T  ) 

E(jUIVAI,ENCF  (Y(??0H)  »GAMAH)  , (Y(P209)  .GAMAE) 

EOUI  valence  (Y  (P0S7)  .LPI  (D) 

EOUIVALENCF (Y (?000) *LRL ( 1 ) ) » (Y (hlO) *VE) 

DIMENSION  IK?) 

REAL  LRL (9) .LPL (9) 

flight  PATH  ANGLES 

In  U( I ) =T 

CALL  ITAH ( 10 . 1 «T .GAMAE ) 

CALL  SENCOS (GAMAH.SH.CH.O) 

CALL  SENCOS (GAMAF .SE»CF.0) 

PPINCTPAL  TO  I OCAI  TRANSEFf?  I»IATRIX 

2n  LPL (l)=CH«rE 
LPL (?) =-SHorF 
LPl (3)=SE 
LPL (4) =SH 
LPL (5) =CH 
LPl (8) =0.0 
LPL (7) =-CH»SK 
LPL  (8)=SH«SE 
LPL (9) =CE 

Al titmof 

30  H=S0RT (XR*»2+YR*»2+ZR“»?) -PE 

LOCAL  VELOCITY  COMPONENTS 

4n  CAI L PATVEr(LRI  ( 1 ) *Y (800) .Y (512) *0) 

VXLE=Y (51?) 

VYLE  = Y (51 3) ♦WF* (RF  ^H) *roS (P51R) 

V7LE  = Y (514) 

VF=SOPT (VXl F»«?*VYLE**?+VZIE**2) 

RE  TURN 
EnO 


TR301280 


TR3014feO 

TR301500 

TR301560 

TR3016in 

TR301b50 


TH30P160 

TR30P170 


TR30?lfl0 

TR302190 

TR302200 

TR302210 

TR302220 

TR30??3() 

TR302?4() 

TR302250 

TR302260 


TR301740 


TR3U1 750 
TR30) 760 
TR301770 

TR302320 

TR302330 
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APPENDIX  G 

MOD3;  3DOF  PARTICLE  TRAJECTORY  WITH  THRUST 


The  purpose  of  this  module  is  to  calculate  the  altitude,  local  velocity 
components,  flight  path  angles,  and  the  princlpal-to-local  transfer  matrix  for 
a vehicle  flying  along  a 3D0F,  point-mass,  ballistic  trajectory. 

The  altitude  is  calculated  from  the  inertial  coordinates  which  have  been 
obtained  through  the  integration  of  the  equations  of  motion; 


Since  the  velocity  of  the  vehicle  with  respect  to  the  inertial  axes  are  also 
known  from  the  integration  of  the  equations  of  motion,  the  local  velocities  with 
respect  to  the  earth's  surface  are  calculated  by  transforming  the  inertial 
velocities  as. 


\e 

'^Y 

II 

Yr 

+ 

o 

+ h)cosi|>j^ 

LE 

. ^LE  , 

^R 

o 

The  total  velocity  with  respect  to  the  earth  is  then. 


The  flight  path  angles,  angles  of  the  velocity  vector  with  respect  to  the 
local  axes,  are  then  calculated  as. 


G-1 


NSWC/WOL  TR  78-59 


The  transfer  matrix  for  transforming  a vector  from  the  principal  axes  to  the 
local  axes  can  then  be  written  in  terms  of  the  flight  path  angles  as. 


= 


COSYjjCOSYg 

-slnYjjCOSYg 

sinYp 


slnY. 


COSY, 


-cosYjjSinYj. 

sinYjjSlnYg 


COSY, 


M0D3  Parameters 


Parameter 


Units 


Location  in  Y Array 


Input 


Output 


feet 

803 

feet 

804 

2r 

feet 

805 

feet 

3000 

^R 

feet/second 

^R 

feet/second 

^R 

feet /second 

“e 

rad/sec 

radian 

h 

feet 

3013 

\e 

feet/second 

512 

feet /second 

513 

LE 

V. 

feet/second 

514 

"LE 

'^E 

feet/second 

610 

degree 

2208 

’'^E 

degree 

2209 

G-2 

2057-2065 

n n n ri  ^ n n n n ri  onnoo 


( 
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SUBROUTINE  M0n3 

POSITION  MOPULE  KOR  A 3D0F  BALI  ISTIC  PARTICl.t  TRAJ. 
3/10/78 

COMMON  Y(4«)40) 

FOUIVAIENCF  (Y(3013)tM) 

EOUIVALENCF ( Y (?0S7)  .LPl  ( 1 ) ) 

equivalence (Y(803).XR).(Y(«04).YR».(Y(P05).7P) 
equivalence  (Y(301S)*RSTR)»(Y(3ni4)  ,TALIR) 
equivalence (Y (3001 ) *WE ). (Y (3000) , RE) 

E(HI  I valence  (Y(S12),VXLE).(Y(S13)«VYIF),(Y(B14) . VZLF) 
EQUIVALENCF ( Y (2?0b) •RAmAH) , (Y(220V) .PAmaE) 

EOUIVALENCF  ( Y (2000  ) «L  RL  ( 1 ) ) • ( Y (f>l  0)  . VE  ) 

REAL  LPL (9) .( RL (R) 

altitude 

m H=S(JRT  (XR*»?»YR**2  + ZR*»?) -RE 

LOCAL  VELOCITY  components 

?0  call  MATVEO (LPL ( 1 ) .Y  («n0)  . Y (SI?)  .0) 

VXLE=Y (512) 

VYLE  = Y (513)  ♦WE»  (RF*H)  «r()S  (PSIR) 

V7LE=y (614) 

VF=SQPT(VX(F**2+VYLE»*?+VZIF»*?) 

FLKiHT  PATH  ANGLES 

30  CAM  ARKTAmcWLF.VXLE  .GAMAH»0) 

77Z=SQRT (VvlF**?+VYLE«*?) 
call  ARKTAN(V7LE .ZZZ*GAMAE .0) 
call  SFNC0S(GAMAH.SHtCu*0) 

CALI  SENCOS (GAMAE ,SE *CE  , U ) 

principal  to  LOCAt  transfer  matrix 

40  LPL(i)=CH«re 
I R(  (?) =-SH«CF 
LPL(3)=SE 
lE'L  (4)=SH 
LPL  (5) =CH 
LPL (6) =0.0 
LPL  (7) =-CH*SF 
I PL (8) =SH«SE 
LPl (9)=CE 

RETURN 

END 


TM301300 

TR301330 


TR30?530 

TH30?570 

TP30P640 

TR302730 


TR30?elO 


TR30?e?0 

TR302830 

TR302840 


TR303100 
TR3031 10 
TH303220 
TR303230 


TH30324O 

TH303250 

TR3U32b0 

TR303270 

TR303280 

TR303290 

TR303300 

TR303310 

TR303320 

TR303380 

TR3033Q0 
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APPENDIX  H 

IMOD4;  6DOF  INITIAL  DIRECTION 
COSINE  MATRIX 


The  purpose  of  this  routine  is  to  calculate  the  initial  direction  cosine 
matrices  for  transforming  vectors  from  the  inertial  to  the  local  to  the 
principal  axes  for  a 6D0F  trajectory  simulation  over  a spherical  rotating 
earth.  The  inertial  to  local  transform  matrix  is  defined  as, 

r cosiJjj^  sin>l)j^sinTj^  -sinii<j^coST^ 


= 


0 

sinil' 


R 


COST, 


-cosi)^  sintp 
R K 


sinx. 


COS4)gCOSTj^ 


and  the  local  to  principal  transform  matrix  as 


f\p^  = 


cose,,cosYw 

M M 

cosi-.sinv.,  - sin(ti,,slne.,cosY,, 
M M M M M 

-simb.-siny,,  - cos(ti,,sine„cosy,, 
M M M M 


M 

cos4.^cosYj^  + sin<t,j^sinej^sinY„  sin^j^ose^ 

-sinCj^cosY^  + cos<^j^slnE^sinY„  ‘^os^.j^cose^ 


The  angles  and  t_  are  the  latitude  and  longitude  angles  in  degrees  and  the 
angles  (I’m  position  angles  of  the  principal  axes  with  respect 

to  the  local  angles.  These  angles  are  shown  in  Figures  5 and  6. 

The  inertial  to  principal  axis  transfer  matrix  can  then  be  expressed  as 


H-1 
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Input 


IMOD4  Parameters 


Parameter 

Units 

Location  in  Y Array 

"r 

degree 

3014 

'I'R 

degree 

3015 

degree 

2066 

degree 

2067 

degree 

2068 

Output 
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C 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


StiHROUTINE  IMflOft 


Mon  PACKAGE  A iiENfcPAL  PURPOSE  6»iOF  PIIInFO  OP 

UPPUinFO  TpaJECTDRV  PRONRA^■ 

3/10/78 


SIXOOISO 

SIXOOlKO 

PIXOOIHO 

SI xoniPo 
sixoopno 
sixon?iu 


TMIS  pOUTImF  CALCiiLATFS  1 HF  IMTIAI  OIhFCTION  COSINE  RATWirfcS  IR| 
Ll P*LPP 


IT  RFOUIRFS  the  EDLLOwINb  INITIAL  CONPITlONS  ON  COPE  j CONTROL 
CARDS 

PsiR=y(301S)  . I AT  ITiiOfc  (OEfi) 

Ta|)R  = y (30U)  . 1 0N6ITU0E  (OFO) 

GAMAM=Y  (20kF,)  . HtADTNli  ANGLE  (OEG)  OF  PRINC  IPAL  AXES 

FPSILM  = Y (?nF7)  . tlF'/ATION  ANGLE  (DFO)  OF  PFINCIPAL  AXES 

PmIM=Y(206h)  , wOLL  ANGLE  OF  PRINCIPAL  AXES 

Common  ycaqao) 

Eon  I valence  (Y(3rH‘i)«PS|R).(Y(3ni4),TAlip) 

EOllIVALENCE(Y(?noE.)  ,GtMAM)  . (Y(?0S7)  •EP'^lLM)  . (Y(POoR)  .PHIM) 
FOUIVAI  ENCE  < V (?00f>)  «LPt  ( 1 ) ) . ( Y (200P)  .1  PP  ( 1 ) ) t ( V ( 20  ?7 ) . ( l P fl  ) ) 
REAL  LWL  (<3)  .1  RP  (9)  .LLP  (P) 

CALL  SENCOS (PSIR.SP.CP.O) 

CALL  SENCOSCTAOR.ST.CT .0) 

LRL  ( 1 ) =CP 
LPL (2) =0. 

LPL(3)=SP 
LRL (4) =SP»ST 
LPl (S) =CT 
LRL  (Fi)  =-CP*ST 
LRL  (7) =-SP*CT 
LPL (fi)=ST 
LRL  <P) =CP*CT 

Call  SENCOS (gamam.sg.co.o) 

CALL  SFNCOSCEPSIl M.SE .oe,0) 

CALL  SENCOsCPHIm.'^P.CP.O) 

Ll P ( 1 ) =CE«CG 

ll  P (?)  =CP«':G-SP*S>- *CG 

Ll  P ( 3)  =-SP»SG-CP»'-E»CG 

Lt  p (4 ) =-CE*SO 

Ll  P (S)  =CP*C(^»SP*SH>SG 

Ll P (ft) =-SP«CG+CP*4E “SO 

I I P(7)=SE 

I I P(«)=SP«rE 

I I P (P) =CP«CF 

call  MATVEr  (I  I P ( ) ) .LRI  ( 1 ) . LHP  ( 1 ) . 2 ) 

return 

END 


.SIX00220 

SIX00230 

SIX00240 

SIX002SO 

sixonafto 

SIX00270 


SIX00330 

S1X00340 

sixonsso 

S1X00370 
SIX0(I3«() 
SIXO()3PO 
SIXU0400 
SIXOOAIO 
SIX004<'0 
SIX00430 
SIX0044() 
SIX004S0 
SIXO04ft0 
SIXl/0470 
SIXUOAHO 
S] X00490 


SlXOOftPO 
SlXUOft  30 
six()nft4o 


THrs  P4(jj 
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A • 


APPENDIX  I 


M0D4;  6DOF  DIRECTION  COSINE  MATRIX 

The  purpose  of  this  module  Is  to  calculate  the  direction  cosine  matrices 
ip  , ipp,  and  for  a 6D0F  simulation  of  a vehicle  flying  over  a rotating 
s^erical  earth. 

First,  the  latitude  and  longitude  of  the  vehicle  are  calculated  as. 


The  inertial  to  local  transfer  matrix  can  then  be  written  as. 


t^RL^  = 


COSIp. 

0 

sini() 


R 


R 


siniJ)j^sinTjj 


COST. 


-costp  sint 
R K 


-sinij^jjCOSTjj 

sinT„ 


COSlf)  COST 
K K 


In  order  to  calculate  the  position  of  the  principal  axes  with  respect  to 
the  inertial  axes,  it  is  necessary  to  define  the  direction  cosine  matrix  as 
follows. 


i 

t 

-> 

II 

j 

3 

-► 

k 

■o 

k 

m a 

r 

1. 

where 


is  a unit  vector  along  the  inertial  axes 
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and 


j 

-► 

i 

-► 

k 


and 


is  a unit 

vector  along  the  principal  axes 

r 

*12 

*13 

‘21 

*22 

*23 

1 

M 

*32 

*33. 

tV’  = 


When  the  above  are  expanded,  the  equations  look  as  follows: 


“ ‘^ll  *'12 


13  '^R 


*'21  ^R  *"22  ^R  "^  *■ 

^ = Si  S + S2  Jr  * 


23  \ 
33  *^R 


If  these  are  now  differentiated  with  respect  to  time  the  following  equations 
result : 


4- 

s 


t 

+ 

5.  , 

T + 

11 

R 

12 

Jr 

13 

21 

^P 

+ 

CM 

CM 

^R 

5 

23 

t 

+ 

^ + 

So 

31 

R 

32 

■^R 

33 

h 

it. 


If  the  following  relationships  are  substituted  into  the  above  equations 

tp  =:  X tp  = r - q 

tp  .:x  Jp  - 1 tp 

itp  =ux  itp  = q Jp  - P Ip 

and  the  individual  components  are  separated  the  following  nine  differential 
equations  are  formed. 


NSWC/WOL  TR  78-59 


*■21  “ P*'31  " ’^*'11 

51  ° - PSl 

52  " ^*"22  " ‘'*^32 

S2  ' P*^32  " ’^*'12 

52  " ‘*'^12  " PS2 

Ss  “ ^*'23  " ‘^*'33 

53  " P*'33  " ^^^13 

S3  " ‘^S3  " ^*"23 


These  nine  equations  can  then  be  integrated  numerically  in  order  to  define 


the  individual  elements  of  the  direction  cosine  matrix 
can  then  be  calculated  as. 


The  other  matrices 


ft 

803 

ft 

804 

z„  ft  805 

R 

p rad/sec  806 


807 

/ 

808 

2009-2017 

RP" 
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SlIHROUTINE  M0D4 
C 

C 3/10/78 

C 

C MOD  PACKAGF  SIXDG.  A OfnEHAL  PURPOSF  6D0F  RllIOFD  OP 

C UNGUTOFO  TPAJFCTOPY  PROGRAM 

C 

C THIS  ROUTInf  calculates  THF  DIRFCTIOM  COSINF  maTRICPS  LRLtLRPf 
C LIP  FOR  A GOOF  SIMULATION  OVER  A ROTATING  SPHFRICAI  EARTH 

C 

C NOTE"  THE  |RPn(q)  MATRIX  ELEMENTS  ARE  DERIVATIVES  AND  MUST  PE 

C lOENTIFIED  ON  CODF  b AND  7 CONTROL  CAPOS 

C 

common  Y(4040) 

EOUI valence (Y (803) .XP) ♦ (Y (H04) .YR) . (Y IHOS) *7R) 

EOUIVALENCf (Y  t203G) »PH1L) « ( Y (8037) , THE  TALI » ( Y ( 8038 ) . PS IL ) 

FOU I VALENCE (Y (3014)  »TA.lR)  . (Y (3015)  .PSIR) 

EOUIVAl ENCF (Y (POOO) .LPL ( 1 ) ) t (Y (8009)  ,l RP( 1)  ) . ( Y (8018) fLHP(M  1 ) ) 
EOUI  valence  (Y  (80f,)  .PP)  , ( Y (807)  .OP)  t (Y  (808)  .RP) 

FOUI ValENCF ( Y (8087) .LLH (1)).(Y(8039).ILH(1)). ( Y (80S7) ,LPL ( 1 ) ) 
RE  AL  I RL (9) .1  RPD (9)  »IHP(9).LLP(9).LI  R(9)«LPL(9) 

DIMENSION  R(9) 

DEG=1 80./3.141598G535H9 
C 

C LATITUOE  AMD  L0N(4ITUDE 

C 

lO  CALI  APKTAm(-YR.7R.TAUP,0) 

CALL  ARKTAm(XR.  (SORT ( ♦ YR**? ) ) . PS  I W . 0 ) 

20  CALL  SENCns(TAUR.STAH.CTAR,U) 

CALI  SENCOS (PSIR.SPSR.CPSR.O) 

C 

C inertial  To  Lf'CAL.  LOCAL  To  PHINCIPAl,  AND  INEHTIAI  TO 

C PRINCIPAL  TPANSKFR  MATRICES 

C 

LRL  ( 1 ) =CPSR 
LRL  (?) =0. 

LRL (3)=SPSP 

LPL (4) =SPSP*STAP 

LRL (5) =CTAR 

LRL (8) =-CPSR*STAR 

LRL (7) =-SPSR*CTAR 

LPL (8) =STAP 

LRL (9) =CPSo*CTAR 

30  LPPD ( 1 ) =HP*LRP ( ?) -OP*l PP (3 ) 

LRPD ( ?) =PP*LRP (3) -RP*LRP ( 1 ) 

LPPD ( 3) =QP*LRP ( 1 ) -PP*LRP (2) 

LRPD (4 ) =HP*1  RP (5) -OP*L  RP (8 ) 

LRPD (S) =PP*LRP (8) -RP*I  RP(4) 

LPPD(8)=QP»1  RP (4 ) -PP*LRP (5) 

LRPD (7) =RP*LRP (8) -0P«l RP (9) 

LRPD(8)=PP*LRP(9)-RP*I  RP(7) 

LPPD  (9)  =(JP*LRP  (7)  -PP*1  PP  (8  ) 

CALL  MATINNMl  RL .p.l  LP) 

40  CALL  MAT VEC (LRP.Ll R«LLP.?) 

C 

C ORIENTATION  OF  PRINCIPAL  AXES  M/R  TO  LOCAL  AXES 

C 

call  ARKTAm (-LLP ( 7) .LLP ( 1 ) . the  TAl  .0) 

IF(THETAL)  43.44.44 

43  THETALsIHETAI *380. 

44  CONTINUE 

CALL  ARKTAM ( -LL P (h) .ILp(S)  .PhU  .0) 


S1XU1260 

SIX0187O 

SIX01290 
SIXU1300 
S1X0I310 
SIX01320 
SIX01330 
SIX01340 
SIX013S0 
SIXOl 380 
SIX01370 
SIX01380 
SIX01390 
SIX01400 
SIX01410 
SIX01480 
SIX01430 
S1X01440 
SIXOI450 
SIX0I480 
SIX01490 
SIXOI500 


SIX01510 

S1X01580 

SIX01530 

SIX01540 


SIXOlbSO 
SIX01560 
S1X01570 
SIX01580 
SIX0159() 
SIX01600 
SIX01610 
SIX01880 
SIX0lb30 
SI X0184O 
SIXOlbSO 
SlXOlbbO 
SIXOlbTO 
SIXOlbSO 
SlX01b90 
S1X01700 
SIXOl  710 
SIX01720 
SIX01730 
S1X01740 


S1X01750 
SlXOWbO 
SIXOl 770 
S1XU17H0 
SIXOl  791) 


l-S 
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IK(KHIL)  46,48t48 

SIX01800 

46 

PHIL=PHIL*160. 

51XU1810 

48 

CONTINUE 

SIX018?0 

CALL  SFNC05(PHIL»SPH.CPH*0) 

51X01830 

PSIL=ASIN(t  LPI4) >*nt6 

51X01840 

50 

CALL  MATINVd LP.B.LPL) 

SIXOlbSO 

RETURN 

5IX01H60 

END 

51X01870 

I 


1-6 
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APPENDIX  J 


MODS;  TARGET  MODULE 


This  module  provides  for  the  locating  of  a target  with  respect  to  the  vehicle. 
The  target  can  be  either  fixed  with  respect  to  the  earth  or  moving  at  a constant 
velocity.  The  initial  position  of  the  target  is  given  by  the  Initial  altitude, 
h^,  latitude,  and  longitude,  The  velocities  are  presented  as  vertical, 

R„,  longitudinal,  V ; and  latitudinal,  V„  . These  relationships  are  shown 

^ T 

in  the  following  sketch: 


TARGET 


The  inertial  coordinates  of  the  target  are  calculated  as, 

^TR 


-R^cosijj^sinT^ 


= R^cos\1;^cost^ 


- 'Wl-  - 
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where 


The  distance  between  the  target  and  the  vehicle  are  expressed  as, 

/ 2 2 2 
DIST  = ^ + AY  + AZ 

where , 

A = target-missile 

AX  = R^siniJ)  - Xj^ 

AY  = -R  cos>|)sinT  - Y 
T K 

AZ  = R„cosi1;cost  - Z_ 

1 K 
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MODS  Parameters 


Parameter 

Units 

Location  In  Y Array 

Input 

hr 

ft 

1 

^T. 

deg 

2 

1 

deg 

3 

'^T 

fps 

4 

T 

fps 

5 

'I' 

Rt 

fps 

6 

R^ 

ft 

3000 

E 

‘"e 

rad/sec 

3001 

sec 

2862 

ft 

803 

ft 

804 

R 

z„ 

ft 

805 

R 

Output 

X „ 

ft 

10 

TR 

Y „ 

ft 

11 

TR 

Z „ 

ft 

12 

TR 

AX 

ft 

13 

AY 

ft 

14 

AZ 

ft 

15 

DIST 

ft 

19 

J-3 


o o o o r>  o 


I 
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suproutine  Mons 

c 

C 3/10/78 

c 

C TARGET  FIXFD  OR  MOVING  AT  A CONSTANT  VELOCITY 

C 

C THE  REOUIRFO  INPUTS  ARE. 

C HTI  = INITIAL  HEIGHT  OF  THE  TARGET  (ET) 

C TAUTI  = INITIAL  LONGITUDE  (OEG) 

C PSITI  = INITIAL  LATITUDE  (DEG) 

C VTTAU  = LONGITUDINAL  VELOCITY  OF  TARGET  (FPS) 

C VTPSI  = latitudinal  VEI  OCITY  of  target  (FPS) 

C RTDOT  = VERTICAL  VELOCITY  OF  TARGET  (FPS) 

C 

COMMON  Y(4Q40) 

equivalence (Y( 1 ) .MTI ) , ( Y (2) ,TAUTI ) , ( Y ( 3 ) « PSI T I ) 
equivalence ( V (4)  ,VTTAU) « ( Y (S)  fVTPSI  ) . ( Y (6) » RTDOT) 
EQUIVALENCE (Y (3000) *PE ) . ( Y ( 300 1 ) , WE ) 
equivalence (Y (PflGfl) .DEI  T) . (Y(299R) .T) 

EQUIVALENCE (Y (2R62) .T1 ) 

equivalence (Y (7) .RT).(y(rt).TAU).(Y(R)»PSI) 
equivalence (Y (10)  ,XTP).(Y(11),YTR).(Y(12).ZTR) 
equivalence  (Y(  13)  .OX)  , (Y(14)  .DY)  . (Y(1E>)  ,DZ) 

F<,/UI  VAI_ENCE  (Y  (H0  3)  .XR)  . (Y(«04)  .YR)  . (Y(H05)  .ZR) 
equivalence (Y ( ig) .01  ST) 


SIX0I910 

SIX0I920 

SIX01940 

SIX019S0 

SIX01960 

SIX01970 

SIX01980 

SIX01990 

SIX02000 

SIX02010 

51X02020 

SIX02030 

SIX02040 

SIX020S0 

51X02060 

51X02070 

SIX02080 

51X02090 

51X02100 

51X02110 

51X02120 

51X02130 

51X02140 

SIX02150 


INERTIAL  CQORUINATES  OF  TARGE! 

DFG=lH0./3. 1415926S358979 
RT=RE+HTI*PTOQT« (T-TI) 

PSI=P5ITI  + ( VTPSI/RT) *OFG* (T-TI  ) 

CALL  5ENC05(P51.5P.CP.n) 

TAU=TAUTI+ ( VTTAU/ (RT*CP) ) «DE6« (T-T I ) 
TAU=TAU+WE* (T-TT ) «DEG 
CALL  5ENC05(T4U.ST.CT.O) 

XTR=RT*SP 

YTR=-RT*CP*5T 

ZTR=RT*CP*rT 

lEEHTTAL  DISTANCES  BE  T'VEEN  TARGET  AND 


51X02160 

51X02170 

51X02180 

51X02190 

51X02200 

51X02210 

51X02220 

51X02230 

51X02240 

5IX022S0 

MISSILE 


DX=XTR-XR  51X02260 
DY=YTP-YR  51X02270 
D7=ZTP-2R  51X02280 
DTST=5QRT (nx*«2+DY**2*nz**2)  51X02290 
RETURN  51X02300 
END  51X02310 
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APPENDIX  K 

MOD6;  PROPORTIONAL  NAVIGATION  SEEKER  MODULE 


This  module  takes  the  relative  locations  of  the  target  and  missile, 
converts  these  linear  displacements  into  angular  relationships  and  utilizing 
the  laws  of  proportional  navigation  calculates  error  signals  suitable  for  a 
control  system.  The  definition  of  proportional  navigation  is  that  the  angular 
rate  of  the  vehicle  should  be  proportional  to  the  rate  of  change  of  the  llne- 
of-sight  angle.  In  order  to  rotate  the  vehicle,  or  more  appropriately  its 
velocity  vector,  it  is  necessary  to  generate  an  acceleration  at  right  angles 
to  the  velocity  vector.  Based  on  the  turning  radius  of  the  vehicle. 


a 


N 


r 

c 


where  the  following  sketch  shows  the  definitions  of  the  parameters. 


K-1 
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Since, 


0)  X r 


a = 


10 


-a  = 


*Y 


the  radius  of  curvature  is 


r = V/t 
c 


The  normal  acceleration  then  becomes 


or  since  the  definition  of  proportional  navigation  is 

Y =K* 


a 


N 


VK  4' 


The  velocity  component  of  the  missile  at  right  angles  to  the  line-of-sight 
can  be  defined  in  two  ways. 


or 

V = R„  4 

X b 


Therefore, 


and 


and 


Rg4  = Vsin+^Qg 


sin  LOS 


R. 


- ■ 


K-2 
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The  relative  locations  between  the  target  and  the  missile  with  respect  to 
‘the  Inertial  axes  were  calculated  in  the  Target  Module.  These  "position  errors' 
can  be  transformed  into  "body"  related  displacements  by. 


1 

>< 

1 

AXr 

II 

AY„ 

P 

R 

where 


MISSILE 


The  angular  relationships  can  then  be  calculated  as 
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It  is  now  possible  to  express  the  required  normal  accelerations  in  terms 
of  the  angular  displacements.  For  a bi-planar  control  system  these  would  be, 


KV  , „ 

®ZP  Rg 


KV  . , 

a^p  = sinX-p 


If  the  vehicle  had  only  planar  lift  and  roll  capability  these  would  be, 

,2 


KV  , n 
^ZP  = R^ 


(pitch) 


sinXp  (roll) 


The  factor 


can  be  thought  of  as  a navigation  parameter,  NV.  This 


is  usually  tailored  to  fit  the  type  of  seeker  available;  i.e.,  it  should  be 


expressed  in  terms  known  to  the  missile  seeker.  For  an  infrared  system.  Re 


or  the  remaining  distance  to  the  target  is  not  known.  In  that  case  and  in 
others  it  may  be  necessary  to  express  the  navigation  ratio,  NV,  as  constant  or 
as  a function  of  time.  It  may  also  be  necessary  to  have  a different  NV  for 
each  of  the  missile  control  functions;  therefore,  the  program  has  been  written  as 


a^p  = NVP  sin0p 
a^p  = NVP  sinX/p 


a.  = NVR  sin4i_ 
L L 


Each  user  will  have  to  program  the  navigation  ratios  to  suit  his  particular 
system. 
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Input 


Output 


MOD6  Parameters 

Parameter 

Units 

Location  in 

AX 

ft 

13 

AY 

ft 

14 

AZ 

ft 

15 

NVP 

2 

f t/sec 

126 

NVY 

r / 2 

f t/sec 

127 

NVR 

2 

ft/sec 

128 

AZP 

r / 2 

f t/sec 

102 

AYP 

. / 2 
f t/sec 

101 

AL 

2 

ft/sec 

100 

rad 

123 

rad 

124 

K 

rad 

125 

K-5 
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SUBROUTINE  M0n6 
3/10/78 

PROPORTIONAL  NAVIGATION 

THIS  PROGRAM  TAKES  THE  TARGET  OISPLACEmENTS*  CALCUI  ATES  THE 
ANGULAR  DISPLACEMENT,  AND  THEN  PROCESSES  IT  IN  ORDER  TO  ARRIVE 
AT  ERROR  SIGNALS  EOR  THE  AUTOPILOT 

COMMON  7(4940) 

EOUIVALENCE ( Y (?009) ,LRP) 

equivalence (Y ( 1?G) ,NVP) , ( Y ( I?7 ) ,NV Y ) , ( Y (1?8 ) ,NVR) 
equivalence  (Y  ( 100)  ,AL)  , (Y(  JOl)  ,AYP)  , (YdOif)  .A/P) 
EQUIVALENCE(Y(1?3) .PHIF ) , ( Y ( 124) .THETAE) , ( Y ( 1 ?S ) , PS TE ) 

REAL  LRP (9) ,NVP,NVY,NVR 

CALL  MATVEr  (LRP.Y  ( 1 3 ) , Y ( 1 (i  0 ) » 0 ) 

CALL  ARKTAm( Y ( 101 ) , Y ( 1 02 ) , PH  I E , 1 ) 

CALL  ARKTAn(Y(102) .Y(lnn) .THETAE.l) 
call  ARKTAN(Y(101 ) ,Y(lnO)  ,PSIE,1) 

NOTE 

PROGRAM  THE  NAVIGATION  RATIOS  TO  SUIT  YOUR  MISSILE 
NVP=Y (126) 

NVY=Y(127) 

NVR=Y (128) 


A7P=NVP*SIN(THETAE) 
AYP=NVY*SIn (PSTE) 

A|  =NVP*SIN (PHIE) 

RE  TURN 
END 


SllOOllO 

S1100130 
SI  100140 
S1100150 
S1100160 
S1100170 

siiooiao 

S1100190 

S1100210 

S1100220 

S1100230 
SI 1002S0 

S1100270 
S11U0280 
S1100290 
S1100310 
S1100320 
S1100330 
S1100340 
S1100350 
SI  100360 
SI  100370 
S1100380 
S1100390 
S1100400 
S1100410 
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APPENDIX  L 

MOD7;  AUTOPILOT/ CONTROL  MODULE 

This  module  is  representative  of  an  autopilot  or  control  system  for  a 6D0F 
simulation.  In  this  elementary  example,  it  is  assumed  that  the  vehicle  has 
a bi-planar  control  system  and  that  the  vehicle  does  not  roll.  This  module 
demonstrates  how  a lead-lag  network  can  be  Incorporated  for  handling  the 
seeker  error  signals,  how  the  missile  heave  and  pitching  motion  can  be 
incorporated,  and  how  the  actuator  dynamics  can  be  included.  This  module 
has  just  the  two  channels,  one  for  pitch  and  the  other  for  yaw.  The  easiest 
way  to  describe  the  system  is  to  refer  to  Figures  L-1  and  L-2.  The  basic 
input  to  this  module  would  consist  of  two  error  signals  received  from  the 
seeker.  A positive  error  signal  in  the  "Z"  channel  arriving  from  the  seeker 
is  calling  for  a correction  in  the  missile  attitude  such  that  would  cause 
the  vehicle  to  be  displaced  in  the  direction  of  the  positive  axis  of  the 
vehicle.  A similar  arrangement  exists  for  the  yaw  channel. 
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I , 

i FIGURE  L-2  Y CONTROL  CHANNEL 
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I 

I 


I 


Input 


Output 


MOD7  Parameters 


Parameter 

Units 

ERROR  Z 

ZKl 

■^1 

sec 

Z^2 

sec 

sec  2 

Zr 

ft/sec 

qp 

rad /sec 

ZK3 

ZK7 

ZLl 

ZL2 

WNz 

ERROR  Y 
YK5 


!3 

sec 

sec 

YK6 

2 

Yr 

ft /sec 

^^4 

rad /sec 

YK8 

YL3 

YL4 

Location  In  Y Array 


102 

300 

312 

313 

301 
814 
807 

302 
306 

308 

309 

332 

333 

101 

304 

314 

315 


813 

808 

303 

307 

310 

311 
334 
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SUBROUTINE  MOOT 

AUTOPILOT  RASEl «ni-PLANAR  CONTROL  DEFLECTIONS 
3/10/78 

NOTE  THAT  THIS  PHOORAM  REOUIRES  COOf  6 ANU  COPE  7 CONTROI 
CAROS  FOR  A,AD,B.Bn*EO7*EOZDtEOZnD.FOY.E0YO.EOVDO 


S11004S0 

S1100A80 

snooAoo 

SllCOSOO 

SU00460 

SllOOSlO 


COMMON  YIAOAO) 

EOUI VAl  ENCF ( Y ( 100) .ERRORX) . ( Y ( lOl ) . ERROR Y) • < Y ( 102) .FRRORZ) 
EOUIVALENCF(Y(300) .ZKl ) . (Y (301 ) tZK?) ,(Y(302) tZK3) 
equivalence (Y (303) .YK4) • (Y(304) tYKb) . (Y(305) .YKf.) 
equivalence (Y (306) . ZK7 ) . ( Y ( 307 ) . YK8 ) . ( Y (308) .ZL  1 ) 
equivalence (Y (30R) .ZL2) t(Y(310)»YL3).(Y(311)tYL4) 
equivalence (Y (312) tTl),(Y(313)fT2)t(Y(314).T3).(Y(315)tT4) 
equivalence (Y(33?) .WNZ) . ( Y (333) »LZ) t ( Y (334) .WNY) , (Y (335) *1  Y) 
equivalence  <Y (316) .AO) . (Y (317) .A) 

EOUIVALENCp(Y(31») tBD) . (Y ( 319) .B) 

equivalence (Y (806) . RP ) , ( Y ( 807 ) . QP ) . ( Y ( 8 0 8 ) . RP ) 
equivalence (Y (81 2) .Xnu) . ( Y (813) .YOU) t ( Y (B14) ,700) 
equivalence (Y (2000) »LR1 ( 1 ) ) . ( Y ( 2009 ) .1 RR ( 1 ) ) 
equivalence (Y ( 359) .THEIM  ♦ (Y(3b0) .RSID) 
equivalence (Y (329) .tOZO()) . ( Y (330) .EOZD) *(Y(331).E07) 
equivalence (V (336) .EOYOO ) , ( Y ( 337 ) tEOYO) , ( Y (338) .EOY) 
EOUIVALENCE(Y(380) .PMR) , (Y(3ei) .OMP) • (Y(3«2) *RMP) 
equivalence (Y (386) t OA ) . ( Y ( 387 ) t DR ) 

EQUIVALENCE(Y(2nS7) »LR|  (1)) 

REAL  LRL (9) ,l RP (9) .LPL (9) f LZ»LY 


S1100S30 
S1100S40 
SllOOSSO 
S1100S60 
S1100S70 
S1I00580 
S1100590 
SI lOObOO 
SllOOblO 
S1100620 
SI  100630 
S1100640 
SllOObSO 
S1100660 
S1100670 
S1100680 
S1100690 
S1100700 
S1100710 


PITCH  CHANNEL 


EI=ZK1*ERR0R7 
AO= (EI-A) /T2 

E0=(T1/T2)*EI*((T2-T1)/T2)«A 
20  THED=7K2*Eo 


S1100720 

S1100730 

S1100740 

S1100750 


YAW  CHANNEI 


EI=YKS*ERRQRV 

BO=(Er-B)/T4 

E0=(T3/T4)*FI* ( (T4-T3) /T4)oH 
40  PSID=YK6*Eo 

MISSILE  ANOIIIAP  pates  in  PRINCIPAL 

CALL  MATVEC (I RP, Y (812) ,Y (383) *0) 

SO  PMP=PP 

qmP=-OP+ZK3*Y (38S) 

R^(P=RP*YK4*Y  (384  ) 

100  Y(323)=0. 

Y (324) = (THFO-OMP) *ZK7 
Y (325) = (PSin-RMP) «YK8 

LIMITS  ON  ACTIIATQR  SICiMALS 


SI  100760 
S1100770 
51100780 
S1100790 

AXES 

S1100800 
SI  100810 
S1100820 
SI  100830 
SI  100840 
S11008SO 
91100860 


IE (Y(324) .QE.ZL?) 
IE(Y(324).I  E.ZLl) 
IE  (Y  (325)  .(;E.YL4) 
IF(Y(325) .LE.YL3) 


Y(324)-ZI  2 
Y (324) =ZL J 
Y (32S) =Yl 4 
Y(325)=Y13 


SI  100870 
SI  100880 
SI  100890 
SI  100900 


o n 
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AC1UAT0H  DYNAMICS 


EOYDD= 

(WNY**?)  * (Y  (325)  -EOY)  -2.0*LY*FOYO**fNY 

S1100910 

EOZDO= 

(WNZ**?)*(Y(324)-F07)-2.0*LZ*EOZn*KNZ 

SU00920 

DA=EOZ 

S1100930 

DB=EOY 

S1100940 

RETURN 

Si  100950 

END 

S1100960 
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APPENDIX  M 

MODS;  6DOF  FORCE  AND  MOMENT  MODULE 

The  purpose  of  this  module  is  to  calculate  the  external  forces  and  moments 
acting  on  a vehicle  in  a 6DOF  simulation.  It  provides  for  nonlinear  aerodynamics, 
winds,  thrust,  and  relatively  small  control  moments. 


The  velocity  of  the  vehicle  with  respect  to  the  earth  is  calculated  by 
transforming  the  inertial  velocities  as  follows. 


Winds  are  then  introduced  in  tabular  form  as  a function  of  altitude.  The  wind 
velocity,  is  tabulated  as  a function  of  altitude  in  Table  Array  No.  3,  and 
the  heading  angle  of  the  wind  is  tabulated  in  Table  Array  No.  4 (see  Figure  8). 
The  velocity  of  the  vehicle  with  respect  to  the  air  is  then  calculated  as. 


If  the  geometric  axes  in  which  the  aerodynamic  data  were  measured  are  skewed 
with  respect  to  the  principal  axes,  the  velocity  can  be  further  transformed 
as 
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and  the  total  velocity  of  the  vehicle  with  respect  to  the  air  is, 


In  order  to  calculate  the  Mach  number  and  dynamic  pressure  it  is  necessary 
to  calculate  the  flow  peroperties  such  as  speed  of  sound  and  density.  These 
values  are  derived  from  the  1969  U.S.  Standard  Atmospheric  Tables. 

The  angles  between  the  body  and  the  flow  vector  are  defined  as: 


The  aerodynamic  forces  and  moments  are  entered  into  the  program  through  a 
set  of  tables.  These  tables  are  as  follows: 
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AERO  COEFF.  TABLE  ARRAY  NO. 


C^(M,  a ) 

12 

Cy(M,  a, 

13 

C^(M,  a. 

14 

C (M,  a 

) 

15 

C^(M,  a. 

^a) 

16 

C,  (M,  a 
P 

) 

17 

Cp  (M,  a 

6 

) 

18 

C (M,  a, 
m 

*A> 

19 

C^(M,  a. 

20 

C (M,  a 
m 

q 

) 

21 

C (M,  a 
"P 

) 

22 

The  forces  along  the  vehicle  geometric  axes  can  then  be  defined  as 


1 

X 

o 

/ 1 

2 \ 

^XG 

c 

— 

THRUST 

0 

1 

N •-< 

o o 

L 

■(2^% 

] *(A)* 

YG 

Sg 

0 

where 


Sc  = "x 


£i_ 


C = C cos4i.  + C sin4i  + C coscj). 

YG  y A y A y^  2V^  A 


^ZG 


pd  . 

^y  2v7 


P A 

The  aerodynamic  moment  coefficients  are  written  as 
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Parameter 


Input 


/!>X,  AY,  AZ 


Output 


MODS  Parameters 


3006,  3007,  3008 


ft 

803 

ft 

804 

ft 

805 

ft 

3000 

ft /sec 

800 

ft/sec 

801 

ft/sec 

802 

rad/sec 

3001 

deg 

3015 

deg 

3002 

deg 

3003 

ft 

3004 

ft2 

3005 

deg 

604 

deg 

386 

deg 

387 

2007-2035 

2000-2008 

806,  807,  808 

577 

Ib/ft^ 

576 

deg 

572 

deg 

573 

deg 

599 

j 

deg 

574  j 

i 

lb 

550  I 

lb 

551  1 

lb 

552  ; 
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Parameter 

Units 

Location  In  Y Array 

\P 

ft-lb 

547 

ft-lb 

548 

ft-lb 

2. 

GP 

ft-lb 

563-571 

fps 

MODS  TABLES 

575 

Table  Array  No. 

Table 

Units 

3 

V„(h) 

fps 

4 

A„(h) 

deg 

12 

C^(M,  a ) 

13 

Cy(M,  a, 

14 

C^(M,  a, 

15 

C (M,  a ) 

16 

Ci(M,  a, 

17 

(M,a) 

P 

18 

(M,  a ) 

6 

19 

C^(M.  a. 

20 

Sl(M.  a,  4>a) 

21 

C (M,  a ) 
m 

q 

22 

CMp(M.  « ) 

23 

THRUST  (t) 

lb 

24 

m (t) 

slug 

25 

s 

I (t) 

XX 

slug/ ft^ 

26 

I (t) 
yy 

slug/ft^ 

27 

zz 

slug/ft^ 

28 

AX(t) 

29 

AX(t) 

30 

AZ(t) 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SIIHwnilTlNE  MODH 

A (iENFPAl  PURPOSE  MJCF  FORCE  AND  MORFM  'lODULE  FOR  A 
MISSIl  E WITH  MOMENT  COmTHOI  IN  TWO  RIANES 

3/10/78 

TARLFS.  KTAH  ( T)  =VW  (H)  . WIND  SPFEf)  (EPS) 

KTAR  (4  ) =AW  (H)  . WIND  A/lMtITH  (IjE(i) 

KTAB( 1?)=CX (M. ALPHA  HAR) 

KTAR ( 1 3) =CV (M. ALPHA  8AR.PHIA) 

ktar ( ] 4 ) =cZ (M. alpha  RAR.PHTA) 

KTAP ( IS ) =CYP («. A|  PHA  RAR) 

KTAB ( 1 h) =rL (P. Al PHA  RAO,PhIA) 

KTAB ( 1 7 ) =rLP (M. ALPHA  RAR) 

KTAP (1 H) =rLn (M« Al  PHA  HAR) 

KTAR  ( IR)  =:rM  (M,  ALPHA  bARtPHiA) 

KT AH (?n ) =rN (P. ALPHA  bAR.PHIA) 

KTAR  (?1  ) =rMr)  (M,  AL  PHA  RAR) 

KTAR  (??)  =:CNP  (M.  AL  PhA  RAR) 

KTAH  (23)  =TMRLIST  (T)  , (IR) 

KTAH  (24  ) =mS  ( T ) . TOTAL  MASR  . ( R|.  U(- ) 

KTAH ( ?S ) =IXX { T ) « (SLUO-FT*«2) 

KTAH  ( 2b)  = I YY  ( T ) . ( SI  U(T-F  T » A ^ ) 

KTAB(27)=IZ7(T) . (SLUfi-ETA*?) 

KTAH(?R)=nXO(T) 

KTAH(2R)=r)YO(T) 

KTAR  ( 3(1 ) =OZG  ( T ) 

n=REFFRtNCF  IT.  (FT) 

A=REFFRENCF  area  (FT«*P) 

0x0  »i  ocatton  of  ('fom.  axes  w/r  to  afpodyn.  Data  axes 
nvo  * 

OZG  ® (NONOTMFNSIONAL.  FT/o) 

PRO  yaw  ANOI  F (DFO)  HeTwN  GEOM.  AXFS  and  PPIN.  AXF'R 
THG  PITCH  


COMMON  Y(404n) 

CnMMON/TAB/7 (SO ) 

EOUI  VAI  ENCF  (Y(S7P).  GR  A '/).(Y(S77).Vm  A CH) 

FOUI  VAL  ENCF  (Y(RnO).XR|)).(Y(Rl)l)  . YRD)  . ( V ( 802)  , /RU) 

EDO  I valence ( Y ( 30  IS)  .PSTR)  » ( Y (3014)  . TADP) 

FODIVAI  ENCE (Y ( TOO) ) »WF  ).  (Y  (3000)  ,RE ).  (V (2027)  .1  LP  ( n ) 

EODl  VAL  ENCF  ( Y ( 20  00)  .LSI  ( 1 ) ) . ( Y ( 30  1 3 ) . E' ) . ( Y ( Sb3  ) » L GP  ( 1 ) ) 

FOUI  valence  ( Y ( 301  r)  .THE  TA(-;)  .(  Y ( 3020  ) .PSIG)  .(  Y ( 3021  ) .phi G) 

E(JU  I valence  (Y(R0  3).XP).(Y(«0  4).YP).(Y(k0S).ZR) 

EOlJ  I valence  (Y(H0b).P),(Y<H(i7).O).(Y(HUH).H) 
FCUIVALENCE(Y(3004).n).(Y(300S).A) 

FOUI  VALENCE  (Y  ( 300b  ) .OXG)  .(  Y ( 3007  ) .dYG)  .(  Y ( 300w  ) .D’/r,) 

Emu  I VAI  ENCF ( Y (S47 I .Ml  P ) . ( Y (S48 ) .MMP)  . (Y (S4R)  .MNP) 

EDIII VAI  ENCF (Y(SS0).FXR).(Y(Sbl).FYR).(Y(SS?).F7R) 

EODI VAL  ENCE (Y(S12).VX|F).(Y(S13).VYlF).(Y(SlA).V7l>' 

EOUI VAL  ENCE  ( Y (SOO)  .XED)  . ( Y (SOI  ) . YFD)  . ( Y { d02)  . 7ED) 

EOUI  valence  (Y  (20S7  ) .1  PI  ( 1 ) ) . ( Y (20JR)  .1  I P ( 1 ) ) . ( Y (20<iR  ) .1  PG(  1 ) ) 

F ('III  VAL  ENCF  ( Y (S72)  . AL  PHA  ) . ( Y (S73)  .RE  T A ) . ( Y (S7«  ) .PH)  A ) 

FOUI  VAI  ENCE  (Y  (S7S)  -VA  ).(  Y (S7b)  .1)0  ).(  Y ( S20  ).  VA  Al  ) 

FOUI valence ( Y (S2 1 ). VAYI  ) . ( v ( S?2 ) . v a 7|  ) 

EoUIVAl.ENCF  ( Y I 34S1  ) .K  ( 1 ) ) 

EOUI  VAI  E NCE  ( Y ( 2R<.v)  . T ) 

FOUI  valence  (Y(S4A).EX(5).(Y(S4b).FYG).(Y(b4b).LZG) 

EOUI valence  (Y (SR 0)  .thrust ) . ( Y (SRR)  .Al  PP) 


RllOl 030 

SI  101040 
SllOlOSO 
SllOIObO 
SU01070 
SI lOI OHO 
SI IO1090 
SI  101 100 
S 1 1 0 I 1 1 (I 
SllOl 120 
SI  101 130 
SllO) 140 
SllOl ISO 
SllOI IbO 
SI  101 170 
SllOl IHO 
SllOl IRO 
SI  101200 
SI  101210 
SI  101220 


SllOlZbO 
Si  1012  70 
SI  10128(1 


S110132(l 
SI  101 330 


S11013SO 
SI  101380 
SI  1013 70 
SllOliHO 
SI 1013R0 
SI  101400 

SI  101420 
S110143(> 
SI  101440 
SI  1014S(: 
SI  101480 
SI  10 1470 
S11014H0 
S11014R0 
SllOlSOO 
SI  101 81 0 
SI 101b20 
S1101S30 
SI  1018411 
SI lOlSbO 
SllOIbbo 
SI  101b 70 
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E0UIVALENCF(V(5?<>)  tCX)  . (Y(532)  «CV>  . (Y(S33)  tCZ) 
E0UIVALENCF(Y(53*)  .CLP)  . (Y(536)  .CM)  . ( Y (‘>37  ) .CmO ) 
E0UIVALENCF(Y(535) .Cl D) . (Y(539) .CN) 

EOUIVALENCF (Y (614) .CYP) . (Y(61b) .CNP) . (Y(616) .Cl  ) 
EOOIVALENCF  (Y  t51«)  . VM ) , I Y (*>19)  . AM) 

EODI VALENCF ( Y (2037) .THFTAL ) . ( Y ( 30 1 2) .MS ) 

EOUIVALENCF  (Y  (S26)  .VAX(;)  . (Y  (527)  . VAYG)  . (Y  (528)  .VAZfi) 
EOUIVALENCF ( Y(2009) .1 WP) . ( Y (558) .PHl AD) . (Y (604) .OAF) 
EOUIVALENCF ( Y (605) .PD2V) . (Y(606) .ODZV) . (Y(607) .RD2W) 
EOUIVALENCF (Y (386) .OA) . (Y(387) .08) 

EOUIVALENCF (Y(608) .CMOA) . (Y(609) .CNOB) 

EOUIVALENCF lY (611 ) .MLG) . (Y (6l2) .MMG) . (Y (613) .MNG) 
EOUIVALENCF (Y (3009) . I XX ) . ( Y ( 30 1 0 ) . I YY ) . ( Y ( 301 1 ) . I Z7 ) 
EOUIVALENCF (Y (610) .VE) 

OIMENSION  K(49) 

OTMENSION  R(9) 

DIMENSION  U(3) 

RFAL  LLR(9) .LPG(9) .LPL(9) 

RFAL  LRL(9) ,LGP(9) .MLP.MMP.MNP.LLP (9) 

RFAL  LRP(9) jMS.IXX.IYY.IZZ 
RFAL  MLG.MM6,MN« 

840=3.141592653589/ 180. 


S11U1580 
S1101590 
S1101600 
S1101610 
51101620 
SI 1UI630 
51101640 
51101650 
51101660 
51101670 
51101680 
51101690 
51101700 
51101710 
51101720 
51101730 
51101740 
51101750 
51101760 
51101770 
51101780 
51101790 


VELOCITIES 

0 CALL  MATVEC (LRL ( 1 ) .Y (800) . Y (512)  .0) 
H=SORT (XR**2*YR**P*ZR**2)-RF 
VXLE=Y (512) 

VYLE=Y (513) .WE* (8F*h) *rnS(RAD*PSlR) 
VZLE=Y (514) 

VF  = SORT (VXl F**?*VYI  E**?.VZLfc**2) 

CALL  ITAB(3.1 ,H,VW) 

CALL  1TAB(4.1.H.AW) 

Akl=AW«RA0 

60  VWXLE=-VW*rOS( AW) 

VWYLE=VW*SIN(AW) 

VWZLE=0.0 

70  VAXL=VXLE-VWXLE 
VAYL=VYLE-VWYLE 
VAZL=VZLE-VW7LF 

CALL  MATVEC (Y (2027) .Y (520) . Y (523) .0) 

PRINCIPAL  AXIS  MISALIGNMENT 

CALL  5ENC05 (THETAG.STG.CTG.O) 

CALL  SENCOS (PSIG.SPG.CPG.O) 

CALL  SENC05(PHIG.SPHG.rPHG,0) 
LGP(1)=CPG*CTG 

LGP ( 2) =CPHG*STG-SPH6*SPG*CTG 
LGP (3) =-SPHG*ST6-CPHG*SPG*CTG 
LGP(4)=-CPG*STG 
LGP (5) =CPHG*CTG.SPHG*SPG*STG 
LGP(6)=-SPG*CTG.CPHG«SPG*STb 
LGP(7)=SPG 
LGP(H)=SPHG*CPG 
LGP(9)=CPHG*rPG 
CALI  MATINVd  GP.P.LPG) 

00  CALL  MATVEC(LPG.Y(523) .Y(526) .0) 

ATMOSPHERIC/FLOW  PROPERTIES 


51101800 
51101810 
S1101820 
51101830 
51101840 
51101850 
51101860 
51101870 
51 loieeo 
51101890 
51101900 
51101910 
51101920 
51101930 
51101940 
51101950 


51102070 

51102080 


51102090 

51102100 

51102110 

51102120 

51102130 
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» 


OP=O.S*HHn»VA**? 

Gn  TO  IIH 
115  VKA(,H  = n.O 

RHO  = 0.() 

OP=0.0 

IIB  GPAV=(32.l74*hlF**?)/((S(JHT  »?>)«*?) 

C 

C ANGULAP  REi  ATIONSHIH  HETwftK  MlSSKt^  AmO  VPLiUITV  wfCIOP 

C 

l?n  CALI  APKTAm  ( Y (S?H)  . Y (S>*5)  . A|  PMA,0  ) 

CALL  APKTAm  ( Y (5?7)  . Y (5.»5)  .hT  TA.d  ) 

Call  abkTam(Y(S?7).Y(s;>h).phia.(i) 

CALI  APKTAm  ( (SOPT  ( Y (S27  ) »*?*Y  **?)).  Y « At  PP.O  ) 

C 

C FORCE  and  moment  GFNFpATIOt!  »tnnni**»*****o*«»*******«*»*«*»»*** 

C 

C TARIILATEO  AFPOOVOAMIC  COKEFICIEmTS 

C 

II  ( 1 ) =\/MACH 
IM?)  =ALPH 

CALL  T TAB  ( 1 ?,?.INrX  ) 
call  TTA8 ( 1 5.P.U.C YP) 

CAl L T TAB ( 1 7 .P.O.ri  R) 

CAI  L TTAH  ( 1 E.R.ll.Cl  U) 
call  T TAB (p)  .P.U.CmO) 

CAt  L T TAR  ( p?.?.Ii,CNR) 

TF(PHTA)  lAO.lSO.lSU 
140  PHIA=PHI A*7hO. 

150  CONTINUE 

U(J)=AMOr)(PHTA.Rn.  ) 

CALI  T TAR ( 1 3. T.U.CY ) 

CALI  ITAH(14.3.(J.C/) 

CALI  1 TAB  ( 1 E>.3.U,ri  ) 
call  TTAB(1R.7,i|,rM) 

CALI  TTAR  (?0.3«INrN) 

C 

C ANGULAR  RATF5  OF  HOL>Y  './R  TO  FLOW 

C 

call  ‘:FNC05(PHTA,sPEt.CPH«t') 

CAl  I MATVEr  (I  PP.Y (Hi?)  .Y (553)  .0) 

CAI  I.  MATVECd  PG.Y(553)  .Y(5G3)  ,0) 
call  MATVEr (I  PR. Y ( HUG)  . Y (H?n ) ,0 ) 

If  (VA7G**2.VAYR«*?)  IGO. 150,170 
150  PHIAnrY(H2n) 

GO  TO  )75 
170  CONTINUE 

PMT  A0=  ( VA7G*  ( Y (554  ) - Y ( J??  ) » v A A'i ) - V A Y(.®  ( Y (555)  , Y ( 5?))  *v  A>G)  ) 

♦ / ( VA7R**?*V  AYG*«?  ) ♦ Y ( »>>>i) ) 

175  CONTINUE 

Pivv=PF(iAr)«n/(?.o«vA) 

On?V  = Y (B21  ) *1)/  ( ?.r,ovA  ) 

Pri?V  = Y (H22)  *0/  (?.n»v  A ) 

c 

C AFRO  rOEFFrriFNTS 

C 

cxG=rx 

CYG=Cv*CPH*C7*5PH.rYP«rJ|W«rPH 
C7G=C7*CPH-CY»SPH-rYP<iPn?V0SHH 
Cl  G = CL*CLP«Pri?V 

CFG  = CMeCPH  + CN*5PH*(;NP«P0?v"SPH*CM(J«0U?u-C5nA«|)A 
rA.G  = CN*CPh-CM«SPH*(  NP‘»P')?V®CPU*f:M(j4pr)?v*CNOR«0P 

c 

C TF'RUST  and  MA55  PPOPF-IIFS 

C 

U ( 1 ) = T 


SI  102140 
SI  102150 
SI  1 OHl  5(/ 
sn  021  70 
SI 1021H0 
SI  1021  VO 


SH02200 
SI  10221 0 
SI  102220 
SI  102230 
S1102240 

51102250 


SI  102270 
SI  102280 
SI 1022V0 
SI  102300 
SI  10231 0 
SI  102320 
SI  l(i?330 
SI  102340 
SI  102350 
SI  102350 
SI  102370 
SI  1 02JHO 
SI  l(l?3<)0 

51102400 
S1102410 
SI  102420 
SI  102430 
SI  102440 


51102450 
SI  102450 
51102470 
SI  102450 
51 1024VO 
SI  102500 
51102510 
51102520 
SI  1025  )o 
511 02540 
51102550 
51102550 
51102570 
SI 1025H0 


SI 1025Mn 
5110250(1 
SI  1025  1 0 

SI  102530 
5 1 1 02540 


Cl  102540 
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CALL  irAHf>3»l«ll*THHUST) 
call  ITAB(?4.1  .I).*-*.) 

CALI  TTAB(?S«1  .II.IXX) 

CALL  ITAB  (?b.  1 .11.  T YY) 

CALL  ITAB(?7.1.IJ.17Z) 

FnPCES  AND  MOMFMS 

FXG=CXG*OP«A*THWMST 

FYG=CYG*QP*A 

F7G=C7G*QP*A 

call  ITAB  (?Fi.  1 .ii.iiXG) 

CA|  I ITAHlPP.l.l'.OYG) 

CALL  !TAB  ( ->0.  I .1J.'17G) 

Ml  6=  (rLG*CYG*n7(3-r7G*l)YG)  *ijp*A*H 
MMG=  (rMG*C7G*DXG)  *OP*A«l) 

MNG=  (rNG-CYG*GXG)  oOP^A*!! 

call  mat  vet  (I  (iP.  Y (G1  I ) .Y  (S47)  .") 

?0P  CfTMTIMlIfc 

CALI  mATVEC (I GP.Y (S44) .Y (bSO) .0) 
call  mATVECIIPI  . Y (SbO  ) . Y (SS()  ) .0) 

CAIL  MATVEr(||P.Y(SbO).Y  (bSti ) . (I ) 

Ph  TIIP»' 

FAin 


';uo?b6o 
SU02h70 
CllOPhHO 
«:i  1026P0 
BllOPTOO 


G1102710 
Gl 10P720 
SI  102730 
SU02740 
SI  102750 
sn027G0 
SUo?7/n 
SI 1027B0 
snO?7PO 
SI  U'PMOO 
SI I02«in 
SI 102P20 
SU0283n 
Cl 102M40 
S1102MSO 
SllOPHbO 


„,S  PWE  U best 
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APPENDIX  N 

MOD9;  EQUATIONS  OF  MOTION  MODULE 


The  purpose  of  this  module  is  to  express  the  equations  of  motion  for  a 6DOF 
simulation  of  a fairly  general  body.  The  force  equations  are  written  for  integration 
in  the  inertial  frame  and  the  moment  equations  for  integration  in  the  principal 
axes. 

i 

For  a mass,  m,  with  the  following  characteristics. 


where  Xp,  Yp,  Z are  the  principal  axes  of  the  mass,  and  X^,  are  the 

inertial  axes,  the  general  acceleration  of  point,  p is. 
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If  the  assumption  is  made  that  m is  a rigid  body,  then 


(ap)  = 0 
(Vp)o  = 0 

and  if  0 is  the  center  of  gravity  of  m,  then  the  following  force  and  moment 
equations  can  be  written. 

F = (a)  dm  - n>(a  ) 

m r OR 

M = ^ [OP  X (a)„]dm  = ^dl*x[uxOP  + (ux(u)x  OP)]dm 
m R m 

These  equations  can  then  be  expanded  into  the  following  forms  for  the  case 
where  the  cross  products  of  inertia  are  zero. 

^R  “ “S  ^XR 

^R  ^ “s  ^YR 
^R  ° “s  ^ZR 


CdLt 


where  and  are  the  resultant  forces  in  the  inertial  system,  and 

M^p,  are  the  resultant  moments  about  the  principal  axes.  In  the  force 

equations  the  resultant  forces  can  be  separated  into  those  due  to  gravity  and 
those  due  to  all  other  forces.  The  force  equations  then  become, 

^XR 

Xr  = ^ - g sin^p 

^YR 

^R  = ^ ^ cosq^psintp 


S - 8 COS.(.pCOSTp, 


where  F F , F „ consist  of  all  external  forces  except  those  caused  by 
, XK  YK 
gravity. 
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Input 


Output 


Parameter 


MOD9  Parameters 

Units 


Location  in  Y Arrav 


degree 

3015 

^R 

degree 

3014 

"XR 

pound 

550 

^YR 

pound 

551 

^ZR 

pound 

552 

•"s 

slug 

3012 

GRAV 

f eet/second^ 

578 

P 

rad /sec 

806 

q 

rad/sec 

807 

r 

rad/sec 

808 

^XX 

2 

slug- ft 

3009 

Iyy 

2 

slug-f t 

3010 

^zz 

2 

slug-f t 

3011 

\p 

feet /pound 

547 

feet/pound 

548 

feet /pound 

549 

^R 

2 

feet/second 

812 

"r 

2 

feet /second 

813 

fR 

2 

feet/second 

2 

814 

P 

rad/sec 

815 

q 

rad/sec 

816 

r 

rad/sec^ 

817 

ODD  non 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUPHOIITINE  MOnQ 


3/10/78 

EOUATTONS  OF  MOTinM 


MOD  PACKAGF  SIXOP,  A (-.FNEPAL  PURPOSE  80UF  GUlOFO  Oo 
unguioed  trajectory  program 

THE  mass. MS  ANO  moment*;  OF  INERTIA  1XX.IYY.I7Z  ARE  PEOUIRFO 
NOTE"  DERIVATIVES  XRliO . YROU » ZRDD . XRD . YRD. 2R0 . PD .UD . RD.R .0 .R 


COMMON  Y(4940 

) 

EOUIVALENCF ( Y 

(S47) 

.MLP)  . ( 

Y (548 

) .MMP)  . ( Y 

EOUIVALENCF (Y 

(SSO) 

.FXR)  , ( 

Y(551 

) .FYR)  . ( Y 

EOUIVALENCF (Y 

(812) 

.XPOO)  . 

(Y(H1 

3) .YRPO)  . 

EOUIVALENCF (Y 

(815) 

.PO)  . (Y 

(818) 

.00) 

. (Y(8 

EOUIVALENCF (Y 

(3nOR 

) • IXX)  . 

( Y (30 

10)  . 

lYY)  . 

EOUIVALENCF (Y 

(808) 

.P) . ( Y ( 

807)  . 

0)  . ( 

Y (m08 

EOUIVALENCF (Y 

{ 301  2 

) .MS) 

EOUIVALENCF ( Y 

(3018 

) .PblP) 

. ( Y (3 

014) 

.TAUP 

EOUIVALENCF (Y 

(2RRR 

) .T) 

RFAL  MS.MLp.MMP.mnP. I XX . I YY. IZ7 


(SaR) «mNP) 

(5S?) .E  ZR) 

(Y (8Ia) .ZROO) 
17) .RD) 

( Y (301 n » IZ7) 

) .R) 

) . ( Y (S781 .GRAV) 


LINEAR  EOUATIONS  OF  MOTION 

CALL  SENCOS(PSIR.SPS.CPS.O) 
CALL  SENCOC  (TAIJR.STA.CTA.O) 
XPOD=FXH/MS-GRAV*SPS 
YPOD=FYR/MS.GRAV*CPS*STA 
ZRDO=FZH/HS-GRAV*rPS«CTA 

ANGULAR  equations  OF  MOTION 

30  PP= (0*R* ( IVY-IZ7) .MLR) /IXX 

00= (P*R* ( 177-TXX ) .MMP) / I YY 
Rn=(Uop*( I*X-IYY) .MNP) /IZZ 
40  CONTINUE 

RF  TURN 
END 


SIX0A4S0 

SIX064f.« 

SIXOGAHO 

SIXOG4RO 

SIXObSOO 

SIXOAblO 

SIX0f.S?0 

SIXOf.530 

SIX06540 

SIX06550 

SlXOftbbO 

SIX0F.S70 

S1X06S80 

SIX0AS90 

SIX0A600 

SIXUFvblO 

SlX0F>6a0 

SIX06630 

SIXOf.640 

SIXOf.f>SO 

SIX0G660 

SIX0F.fr70 

SIXOGhHO 


SIXOGGRO 

SIX06700 

SIX0A710 

SIX06770 

SIX06730 


SIXOG740 

SIX0f>7St) 

SIX06760 

S1X0G770 

SIX0h780 

SIX0G790 
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APPENDIX  O 

MOD14;  3D0F  FORCE  AND  EQUATIONS 
OF  MOTION  MODULE 

The  purpose  of  this  module  is  to  calculate  the  forces  acting  on  the  vehicle 
flying  along  a particle  (3D0F)  trajectory  and  to  set  up  the  equations  of  motion. 

If  winds  are  desired,  they  may  be  Input  in  tabular  form.  The  wind  velocity, 

V , and  the  heading  angle,  A (measured  clockwise  from  the  north)  can  be  tabulated 
as  a function  of  altitude,  h'^in  TABLE  ARRAYS  NO.  3 and  4 respectively.  In  that 
case  the  velocity  of  the  vehicle  with  respect  to  the  air  is. 


1 

< 

f 

1 

w 

1 

-V  cosA 
w w 

\l 

= 

LE 

- 

V sinA 
w w 

0 

_ — 

— — 

or 

V =-/v  ^+V  ^+V  ^ 

^ y ^L  A^L 

The  atmospheric  properties  such  as  the  speed  of  sound  and  density  are  gotten 
from  the  1969  Standard  Atmospheric  Prrpertles  Tables. 

The  forces  along  the  body  principal  axes  are  then  expressed  as 


C is  the  drag  coefficient  input  in  tabular  form  as  a function  of  Mach  number  in 
TABLE  ARRAY  No.  5. 


0-1 
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Other  associated  necessary  quantities  are  the  mass  of  the  vehicle  and  the 
acceleration  due  to  gravity.  These  are, 

nig  = nij  + m(t  - t^) 

where 

iDj  = Initial  mass 
m = mass  depletion  rate 
t^  = Initial  time 
mg  = system  mass 
t = current  time 

and 


The  equations  of  motion  can  then  be  written  as. 
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Input 


Output 


MOD14  Parameters 


Parameter 

Units 

Location  in  Y Array 

H 

« 

ft 

3013 

fps 

512 

LE 

"'y 

LE 

fps 

513 

V. 

fps 

514 

^LE 

ft 

803 

Y„ 

ft 

804 

R 

Z„ 

ft 

805 

R 

slug 

658 

m 

slug/sec 

589 

sec 

2863 

t 

sec 

2999 

THRUST 

lb 

590 

V. 

fps 

520 

hh 

fps 

521 

\l 

% 

fps 

522 

ZL 

fps 

575 

M 

fps 

577 

s 

fps 

529 

F 

lb 

601 

XP 

P 

lb 

602 

YP 

F _ 

lb 

603 

ZP 

F 

lb 

550 

XR 

F 

lb 

551 

YR 

F 

lb 

552 

..ZR  0 

ft/sec 

812 

..R 

Y„ 

2 

ft/sec 

813 

..R 

Zp 

2 

ft/sec 

814 

0-3 


r>nn  rjoo  n rt  n n ooonn 
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SUBROUTINE  MOO 14 
C 

C FORCE  MOOULE  FOR  A 300F  PARTICAl  TRAJ. 

C 

C 3/10/78 

C 


c 

TABLES.  KTAB(3)=VW(H)  . WIND 

SPEED  (FPS) 

TR301 360 

c 

KTAP (4) SAW (H)  . WIND 

azimuth  (DEG) 

TB301370 

c 

KTAP(5)sCn(M)  .DRAG 

COEFF. 

c 

c 

KTAR(23)sTHRUST(T) 

(LR) 

c 

c 

INPUT  parameters 

C 

C 


Mr  = Y(f,58)  , 
MnOT  = Y (589)  , 
TT=Y(2862) 
THRUST=Y(590) 
A=Y(3005)  , 


INITIAL  MASS  (SLUG) 

RATF  OF  CHANGE  OF  MASS  (SLUGS/SFC) 
initial  time  (SEC) 

• CONSTANT  THRUST  (LH) 

RFF.  AREA  (FT**2) 


COMMON  Y(4Q40) 


EOUIVALENCF 

(Y 

(3013 

) .H)  . 

EOUIVALENCF 

(Y 

(512) 

.VXLF 

EOUIVALENCF 

(Y 

(2057 

) .LPL 

EOUI VAl ENCF 

(Y 

(803) 

.XR)  . 

EOUIVALENCF 

( Y 

(578) 

.GRA  V 

E(5UIVALENCF 

(Y 

(812) 

.XPDO 

EOUIVALENCF 

(Y 

(3012 

) .MS) 

EOUIVALENCF 

(Y 

(550) 

.FXR) 

EOUIVALENCF 

(Y 

(3005 

) .A)  . 

EOUIVALENCF 

( Y 

(349) 

.PSIR 

EOUIVALENCF 

(Y 

(575) 

.VA)  . 

EOUIVALENCF 

(Y 

(521  ) 

.VAYI 

EOUIVALENCF 

(Y 

(601  ) 

.FXP) 

EOUIVALENCF 

(Y 

(589) 

.MOOT 

EOUIVALENCF 

( Y 

(529) 

.CD) 

EOUIVALENCF 

(Y 

(658) 

.MI ) . 

DIMENSION  II 

(2 

) 

REAL  LLR(9) 

.L 

PL  (9) 

. MS  . M 

Y(2999) tT) 

» (Y(513)  .VYLF)  f (Y(S14)  .VZI.F) 
1)  ) * (Y(2039)  .LLRd)  ) 

Y (804) ,YR) t ( Y (805) .ZR) 

» ( Y (577) « VMACH) 

* ( Y (813) . YROD) . ( Y (814) .ZROO) 

(Y(551) tFYR) , (Y(552) .FZR) 
Y(3000) tRE) 

♦ (Y (348) *TAUPI 
Y (576) .OP) t ( Y (520) .VAXL) 

. (Y(522)  .VAZL) 

(Y(602) .FYPJ . (Y(603) .FZP) 

» (Y (590) .THRUST) 

Y(2862) .TI) 

.MOOT 


TR301450 

TR301480 

TR301520 

TR301530 


TR301580 

TR301590 

TR301620 

TR301630 

TR3U1640 

TR301670 


PAO=3.1415Q2653589/1BO. 


TR301820 


WINDS 


call  TTAB (3. 1 .H.VW) 
CALL  ITAB (4. 1 .H.AW) 
AW=AW»RA0 

60  VWXLE=-VW«rOS(AW) 
VWYLE=VW*STN(AW) 
VWZLE=0.0 


TR301830 

TR30I840 

TR30ie50 

TR301860 


VELOCITY  W/H  TO  AIR 


70  VAXL=VXLE-VWXLE  TR301870 

VAYL=VYLE-VWYLF  TR301880 

VAZL  = VZLE-VWZI.E  TR30ie90 

VA=S(3RT ( Y (520) **2*Y (521 ) **2.Y (522) **2)  TR301900 

ENVIRONMENTAI  PROPERTIES 

IF (H.GE. 500000. ) GO  TO  115  TR301910 

110  CALL  ARDCFT (H.PP.TT.DO. VS.G)  TR301920 

0-4 


ooo  ononno 
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VMACH=VA/ (VS*1 1 16.4) 

RHO=0n*0. 0023769 

TB301930 

OP=0.5*RHO*VA**2 

TR301970 

GO  TO  118 

TR301980 

115 

VMACH=0.0 

TR301990 

RHO=0.0 

TR302000 

op=o.n 

TR302010 

llfl 

C 

C 

GPAV= (32. 174*RE**?)/ ( (SORT ( XR**2* YR**2*ZR**2 ) ) **?) 

MASS 

TR302020 

TR302030 

TR30P080 


FORCE*; 


U(1)=VMACH 

CALL  ITAB(S.1.U*CD) 

FXP=THRUST-QP*A*CO 

FYP=0.0 

F7P=0.0 

CALL  MATVEr (LPL . Y (601 ) . Y (5S0) *0 ) 
CALL  MATVEr(LI  R.Y(550) ,Y(550) .0) 


FOUATTONS  OF  MOTION 


XRnO=FXR/M«;-GRAV* 

YROO=FYR/MS+GRAV* 

ZROn=F7R/MS-GRAV« 

RETURN 

END 


(SIN (PSIR) ) 

(COS (PSIR) *STN (TAUM) ) 
(COS(PSIR>*C0S(TAUH) > 


TR302090 

TR30?130 

TR302U0 

TR302150 

TR30P270 

TR3OP20O 


TR302290 

TR30P300 

TR302310 


I 
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pr(x:ess  parameters 


Input 

Parcuneter 

Units 

RE 

FT 

WE 

rad/sec 

TAUR 

deg. 

PSIR 

deg. 

TI 

sec . 

T 

sec . 

VXLE 

fps 

VYLE 

fps 

VZLE 

fps 

Output 

GX 

g 's 

GY 

g's 

GZ 

g's 

GAMAH 

deg. 

GAMAE 

deg. 

Location  in  Y Array 

3000 

3001 

3014 

3015 
2862 
2999 

512 

513 

514 


2205 

2206 

2207 

2208 
2209 


non  r>  n r>  r>  n r>  non  n n n n r>  n 
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SUBROUTINE  PROCESS 
3/13/70 

MOD  PACKAGE  SIXDG*  A GENERAL  PURPOSE  GDOF  GUIDED  0® 

UNGUIDED  TPAJECTORV  PROGRAM 

COMMON  VIAQAO) 

EOIII VALENCE (Y (?863 1 t STOP  I » (Y (3012) .MS) t ( Y(6S«) tMI ) 
E0UIVALENCE(Y(2213) .OT ) . < Y (221 4 ) tOP) t < Y ( 221 5 ) .RS ) 
equivalence (Y (3000) »RK) . ( Y (20b2) f T I ) « ( Y ( 300 1 ) , WE ) 

EQUIVALENCE (Y (504) .01 Al. ) . ( Y (505) .OLP) » (Y(506) .OTIM) . (Y(299Q) . I ) 
EQUIVALENCE (Y (3014) .TAUR) . ( Y ( 30 15 ) .PSl R ) 
equivalence (Y (81 «) .RTE I ) » (Y (019) .RPEI ) 

EQUIVALENCE (Y(803) .XRM) . (Y(804) .YHM) • (Y(805) .2RM) 

EQUIVAI.ENCF  (Y  (2200)  .PTE  ) « (7(2201  ) .RPf  ) 

EQU I VALENCE (Y (2205) .GX) . (Y (2208) , GY) . ( Y (2207) .GZ) 

EQUIVALENCE (Y (51?) .VXLF) . (7(513). VYLE).(Y(514).V7lE) 

EOUI valence (Y (220M) .GAM AH) . (7(2209) .GAMAE) 
equivalence  (Y(??04).GT) 

equivalence (Y (GSl ) .CXA) . (7(852) .CY A ) . ( Y ( 653 ) .CZ A ) 

EQUIVALENCE (Y (528) .VXG) . (Y(527) .VYG) . (Y(528) .yZG) 

EQUIVALENCE (Y (648) .CXG) . ( Y (649) .CYG) . ( Y (650) .CZG) 

EQUIVALENCE (Y (880) .L6A ( 1 ) ) . ( Y ( 3005) . A ) . ( V ( 578 ) .GRAW) 
equivalence (Y( 3012) .MS) . (7(854) .BALC) 

REAL  MI.MS.LGA(9) 

EQUATORIAL  PANGES 

RAD=3. 14 1502853509/180. 

RTE=(RE*(TAUP*RAn-WE*(T-TI) ) ) 

RPE= (PE*PSTP*PAD) 

BODY  accelerations  IN  THE  PRINCIPAL  AXES.  (G.S) 


CALL  MAT VEC (Y (2000) .7(812) .7(2202) .0) 

6X=Y(2202) /32.174 
GY=Y(2203)/32.174 
GZ=Y(2204) /32.174 
GT=S0RT (GY**2*GZ**?) 

LOCAL  flight  path  ANGLES. (OEG) 

CALL  APKTAm(-VYLE.VXLE.GAMAH.O) 

77Z=SORT ( VxlE**?*vYLE**?) 

CALL  ARKTAN(VZLE.ZZZ.lSAMAt.O) 

true  distance  traveled  w/p  to  THE  earths  SUREaCE 

TAUR=TAUR«pAD 

PSIR=PSIR*PAn 

TAUE=TAUR-WE*T 

DX2= (RE*SIN(0P) -RE*SIN (PSIP) ) **2 

072= (-Rt*CnS(OP)*STN (OT) ♦RE*COS (PS  I P ) *S I N ( T AUF ) ) **? 
D72= (RE*C0S (OP) *COS(OT) -PE*COS(PSIR) *COS (TAUE) ) **? 
DC=SOPT (DX?*DY2*D72) 

ThETA=2.0*ASTN(OC/(2.0*RE) ) 

0PS=RE*THETA 

RS=RS+0HS 

0T=TAUF 

OP=PSTP 


P-3 


SIX08820 

SIX08B30 

$1X08050 

SIXO806O 

SIXU887U 

SIXO808O 

SIX08R9U 

SIXU8900 

S1X08910 

S1X08920 

51X08940 

SIX0h95U 

51X06980 

51X06970 

51X08900 

51X08990 

51X07000 


51X07020 

51X07030 

SIX07040 

51X07050 

SIX07060 

51X07070 

51X07080 

51X07090 

SIX07100 

51X07110 

51X07120 

51X07130 

51X07140 

51X07150 

SIX07180 

51X07170 

51X07180 

51X07190 

51X07200 

5IX0721P 

51X07220 

SIX07230 

51X07240 

51X07250 

SIX07260 

51X07270 

51X07280 

51X07290 

51X07300 

51X07310 
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TAUR»TAUR/RAn 
PSIRsPSIR/PAD 
IF(STOP)  ?o,?o*in 
in  TIrT 
MT=MS 

20  CONTINUE 
RETURN 
END 


SIX07320 

^1X07330 

S1X073A0  > 

S1X07350 

SIX0736C 

S1X07370  • 

SIX073R0 

SIX07390 


I 

F 
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